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Comments by the Editors 


Raymond C. Montgomery 
NASA Langley Research Center 
Hampton, Virginia 

and 

Howard Kaufman 
Rensselaer Polytechnic Institute 
TVoy, New York 

Hie papers and abstracts contained in this report represent both formal pre- 
sentations and experimental demonstrations at the Workshop on Selected Topics 
in Robotics for Space Exploration which took place at NASA Langley Research 
Center, 17-18 March 1993. This workshop grew from discussions between Dr. 
Ray Montgomery from NASA Langley Research Center (LaRC) and Dr. Howard 
Kaufman from Rensselaer Polytechnic Institute (RPI). Because both the Guid- 
ance, Navigation, and Control Technical Committee (GNCTC) of the LaRC and 
the Center for Intelligent Robotic Systems for Space Exploration (CIRSSE) at 
RPI, shared common research directions in robotics for space exploration, it was 
evident that a forum for technical exchange would be very valuable. 

Thus, with approval from both CIRSSE and the GNCTC, Drs. Montgomery 
and Kaufman solicited papers and/or demonstrations from LaRC, CIRSSE, and 
from persons from industry, government, and other universities with close ties to 
either LaRC or to CIRSSE. 

The presentations were very broad in scope with attention given to space 
assembly, space exploration, flexible structure control, and telerobotics. 
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The organizers would like to thank all those who contributed to the success 
of the workshop. Special thanks are due Betty Lawson of CIRSSE, who retyped 
many of the abstracts and who helped with the initial solicitation of papers and 
Mr. Jack Pennington and Dr. Robert Williams who organized and participated in 
the LaRC Automation Technology Research Branch tours. 
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Introduction to CIRSSE 


Alan A. Desrochers, CIRSSE Director 
Rensselaer Polytechnic Institute 
Troy, New York 

Robotic systems combine mechanical devices such as robot manipulators 
and hands; sensors, such as computer vision, optical, t act i le , and force with 
computer hardware and software to form an integrated system capable of reliably 
executing a variety of different tasks. Much of the research in robotics has 
focused on the theory, algorithms, and development of these individual component 
technologies. In the CIRSSE research program, we emphasize the integration 
of these technologies through the definition of architectural principles and the 
development of an integrated testbed for experimental studies and demonstration 
tasks. We have demonstrated assembly of space-truss structures, and have worked 
closely with NASA laboratories in the definition of space mission scenarios and 
applications. 

The CIRSSE integrated testbed includes two 9-degree -of-freedom robotic ma- 
nipulators mounted on a 12 foot long track. The resulting useful work volume of 
the 1 8-degree-of-fr eedom robotic system is more than 400 cubic feet. The effec- 
tive control of cooperating manipulators for coordinated motion and force control 
tasks is an important research area that utilizes both simulation and experimental 
studies to verify new concepts. 

The CIRSSE testbed incorporates several different types of sensors. There are 
currently five cameras integrated into the system. These TV cameras are supported 
by a special purpose high-speed computer for image processing and automated 
inteipretation of images from a single camera or from two stereo pairs of cameras: 
one fixed and one mobile. A laser scanning range sensor is mounted from the 
ceiling of the testbed in order to provide 3-D depth information from any point 
in the robot work space. Tactile and force sensors are mounted on each of the 
robot hands, and they provide sensory feedback for the dexterous manipulation 
of objects and devices. The interpretation of sensory information acquired from 
multiple sensors is an important research topic in integration. 



The CIRSSE hardware computing environment is based on a distributed multi- 
processor system. Two VME-bus systems support the distributed multiprocessor 
environment and partition the two most critical real time functions: motion con- 
trol and vision processing. The host and development environment is based on 
SUN workstations. 

The CIRSSE Testbed Operating System (CTOS) is built on the commercial 
VXworks software, and provides a versatile multiprocessor real-time operating 
system capability. CTOS supports a versatile message passing protocol that adapts 
to different hardware configurations and available communication speeds. The 
communication protocols for intraboard, intrabus, and network co mm unications 
among processors provides for an efficient implementation and deb ugg ing of 
systems. The message passing protocol that underlies all the CIRSSE testbed 
implementations supports an event driven, object-oriented approach to the overall 
architecture. 

The CIRSSE Testbed integrated architecture is structured hierarchically in intel- 
ligence and integrates planning, coordination, and execution functions. Within 
this broad hierarchical framework the event driven, distributed nature of the sys- 
tem is maintained. Based on the task decomposition, the subplanning modules 
are available both on-line and off-line. The task decomposition itself is repre- 
sented as a Petri net, and this Petri net structure is used as an embedded model 
for the coordination of individual execution modules. The Petri net itself may 
be thought of as a distributed event driven system, and the coordination level is 
implemented using a message passing protocol which parallels the lower level 
execution. The architecture at all levels remains modular and therefore, provides 
flexibility, reliability, and ease of implementation. 

The CIRSSE multi-manipulator, multi-sensor integrated testbed, the CTOS oper- 
ating system, and the CIRSSE integrated testbed architecture are unique develop- 
ments of the CIRSSE Center. They provide a broad resource for the expe rimen tal 
study of robotic concepts and the demonstration of application tasks meeting the 
special needs of space exploration. 
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A Global Approach to Kinematic Path Planning to 
Robots with Holonomic and Nonholonomic 

Constraints 

Adam Divelbiss Sanjeev Seereeram John T. Wen 
Department of Electrical, Computer and Systems Engineering 
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Troy, NY 12180 

divelbis@cat.rpi.edu seereeram@ral.rpi.edu wen@ral.rpi.edu 


Abstract 

Robots in applications may be subject to holonomic or nonholonomic constraints. 
Examples of holonomic constraints include a manipulator constrained through the con- 
tact with the environment , e.g., inserting a part , turning a crank , etc. f and multi- 
ple manipulators constrained through a common payload . Examples of nonholonomic 
constraints include no-slip constraints on mobile robot wheels , local normal rotation 
constraints for soft finger and rolling contacts in grasping t and conservation of angular 
momentum of in-orbit space robots. The above examples all involve equality constraints; 
in applications , there are usually additional inequality constraints such as robot joint 
limits , self collision and environment collision avoidance constraints , steering angle 
constraints in mobile robots, etc. 

This paper addresses the problem of finding a kinematically feasible path that sat- 
isfies a given set of holonomic and nonholonomic constraints, of both equality and 
inequality types. The path planning problem is first posed as a finite time nonlin- 
ear control problem. This problem is subsequently transformed to a static root finding 
problem in an augmented space which can then be iteratively solved. The algorithm has 
shown promising results in planning feasible paths for redundant arms satisfying Carte- 
sian path following and goal endpoint specifications, and mobile vehicles with multiple 
trailers. In contrast to local approaches , this algorithm is less prone to problems such 
as singularities and local minima. 

1 Introduction 

Controlling robot motion, including both in a.iii pnla.tnrs and mobile vehicles, usually 
involves the following steps: 

1. Kinematic path planning: find a path that satisfies all the geometric specifications 
and constraints of a given task. 

2. Trajectory generation: index the path with time to generate a dynamic trajectory. 

3. Dynamic trajectory following: design a servo controller (possibly incorporating 
dynamic information of the overall system) to follow the trajectory. 
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This paper focuses on the kinematic path planning problem. In contrast to most of 
the existing algorithms which are local and reactive in nature, our approach is a global 
one which warps an entire path to satisfy all the constraints. A common classification of 

involves 1 hol ° nom ^ versus nonholonomic. If a constraint can be expressed in 
terms of the generahzed coordinate, it is holonomic; if a constraint involves the gener- 
alized velocity and it is not mtegrable, then the constraint is nonholonomic. Examples 
of holonomic constraints include a manipulator constrained through its contact with 

'T TtiRg a Part ’ tUrning a Crank ’ etc ‘’ and manipula- 

tosrs constrained tirough a common payload. Examples of nonholonomic constraints 

mdMe m-Ap on mobile vehicle wheels, local normal rotation constraints 

for soft finger contacts m grasping, and conservation of angular momentum of space 
mW There may be other design constraints imposed by the task, for example 
equality constraints such as the desired terminal configuration, specified end effector 
path, etc., and inequality constraints such as manipulator joint limits, self collision and 
environment collision avoidance constraints. 

There is abundant literature on path planning for redundant robots, which are 
examples of systems with holonomic constraints, and mobile robots which are examples 
of systems with nonholonomic constraints, but seldom on both. The principal reason 
is that in the local approach, the two problems are fundamentally different in that 
redundant robots can in general move in all directions locally in the configuration 
space, but nonholonomic systems can only move in certain directions. Consequently 
the issues related to redundant robots are singularity, redundancy resolution, joint 
cyclicity for cyclic end effector paths, etc., and, for nonholonomic systems, the main 
issue is in finding a path whose tangent lies within the admissible directions. There 
are also commonalities in the two classes of problems: 


1. Kinematic models are linear in control (i.e., admissible velocities). 

2. Collision avoidance and joint limits are represented as a set of inequality con 
straints. 


In our approach, the entire path is iterated toward a solution. Therefore, in this 
framework, whether the constraint is holonomic or nonholonomic does not make a 
fundamental difference; a more important distinction is between equality and inequality 
constraints. Indeed, path planning for redundant robots and mobile robots are treated 
in exactly the same way in our proposed approach. 

The manipulator path planning problem is traditionally based on geometric ap- 
proaches. (Comprehensive surveys can be found in [1, 2, 3].) The majority of these 
use the configuration representation of the manipulator and the obstacles and joint 
limits (as originally proposed in [4]). In this formulation, the path planning problem 
is reduced to finding a feasible path for a single point. Various tools from geometry, 
topology and algebra have been utilized to develop methods such as roadmap, cell 
decomposition, and potential field methods. The drawback of the configuration space 
approach is that the configuration space can become high dimensional for manipu- 
lators with several joints and mapping out the free region for a redundant arm in a 
cluttered environment is a time consuming process and produces large graphs to be 
searched. Recently, some work has emerged which focuses on minimizing this growth 
in computation by incorporating heuristics [5, 6]. 
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The specific issue of redundancy resolution in redundant manipulators has tradi- 
tionally been addressed by using the Jacobian pseudo-inverse [7, 8, 9, 10], which is 
a local, or point-wise on the path, approach. This approach is simple to use, has 
low computational overhead, and can be combined with the local potential field for 
collision avoidance. The drawbacks include local minima, breakdown near or at sin- 
gularities, and non-cyclic or even unstable joint motion for cyclic end effector paths 
[11, 12, 13, 14]. Some global techniques based on optimization have also been proposed 
[15, 16, 17, 18, 19], but the resulting two-point boundary value problem (TPBVP) is 
difficult to solve for complex problems. 

In the rapidly accumulating literature on nonholonomic motion planning, there 
appears to be a similar dichotomy of approaches: graph search based methods as 
in [20, 21, 22, 23, 24, 25, 26] and analytic methods [27, 28, 29, 30, 31, 32, 33, 34, 
35]. The former class tends to be computationally inefficient but can handle general 
constraints, the latter class gives elegant insights into the structure of the solution but 
may generate impractical paths. Recently, some important insights have been gained 
by linearizing the kinematics equation about a non-stationary trajectory instead of 
a fixed equilibrium. It is shown in [36] that the linearized time varying system is 
frequently controllable, and a locally stable feedback controller can be designed. In 
[37], a globally stable feedback controller is also found. Based on this approach, a 
general procedure for designing globally stable time varying feedback controller has 
been obtained [38], 

A common starting point for the analytic approaches to the path planning problem 
is to pose it as a general nonlinear control problem. Using this framework, we present 
a new method of solving the path planning problem in this paper. The formulation is 
based on converting the nonlinear control problem into a nonlinear algebraic equation. 
The problem then becomes a nonlinear root finding problem in which the dimension 
of the search space is very high (infinite if (u(t) : t £ [0, 1]} is taken from an infinite 
dimensional functional space) compared to the number of equality constraints (for 
the end point constraint). The nonlinear root finding problem is further converted to 
an initial value problem (IVP) which is much easier to solve than the usual TPBVP 
typically arising in the optimization approach. Under some additional assumptions, 
the IVP can be shown to be wellposed [39] and a variable step size ODE solver is used 
to propagate the solution. Examples ranging from a front-wheel driven car to triple 
trailers to nine degree-of-freedom (DOF) manipulators have been successfully tackled 
[40, 41, 42]. A similar approach has also been proposed independently in [43, 44] 
for kinematic path planning with only the end point constraint. For this case, the 
wellposedness property of the IVP is shown to be generic. 

Since the dimension of the search space is very high, there are many possible solu- 
tions to the root finding problem which results from the path planning problem. For 
any practical applications, additional constraints must be placed. We have adopted an 
approach similar to the global exterior penalty function method [ 45 , 46 ] which converts 
inequality constraints into a zero finding problem. This differs from the familiar arti- 
ficial potential field method which is an interior penalty function (or barrier function) 
method in that the initial guess may be infeasible. 

The inequality constrained case cam then be combined with the equality constraints 
into an augmented zero finding problem. Non-configuration space constraints, for 
example, constraints involving corners of a vehicle, body of the robot, etc., cam also be 
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incorporated in this formulation. 

The global path p lanning approach that we have proposed has been applied to 
examples involving redundant manipulators and nonholonomic vehicles. We note the 
following attractive features of this algorithm: 

1. Both equality and inequality constraints can be included in this formulation. 
The constraints can be nonlinear in the configuration variable, so task space 
constraints (which involve nonlinear kinematic function of the configuration vari- 
ables) are also allowed. 

2. The initial guess does not have to be feasible. The planner iteratively warps the 
path until all constraints are satisfied. 

3. This approach emphasizes feasibility over optimality in contrast to other global 
approaches. Once a feasible solution is found, optimality can be incorporated as 
a secondary constraint. 

4. The IVP formulation is computationally easier to solve than the TPBVP. 

5. There are additional points related specifically to redundant manipulators: 

• Goal task variability: This approach can be used for finding a feasible joint 
sequence from a fixed initial configuration to a specified Cartesian or joint 
space goal - the path planning problem - as well as for global redundancy 
resolution along a specified Cartesian path. 

• Singularity robustness: The global nature of the planner avoids the Jaco- 
bian singularity problem inherent in local methods. While the controllabil- 
ity about a configuration is lost at the singularity, the algorithm can proceed 
as long as controllability about the planned path, which is a much looser 
condition to satisfy, is retained. 

• Incorporation of additional equality constraints: As an application of this 
capability, the planner can generate cyclic joint space motion for a specified 
cyclic task space motion. 

There are also additional points related specifically to nonholonomic systems: 

• The planner only requires local controllability about a path in every iteration 
but does not require controllability about a configuration. This is important 
since the latter is not satisfied for nonholonomic systems. 

The rest of this paper is organized as follows: Section 2 describes the theory be- 
hind the proposed algorithm. Section 3 describes the global exterior penalty function 
method to handle inequality constraints. Section 4 shows a number of simulation ex- 
amples involving path planning for redundant manipulators and mobile vehicles with 
trailers. 


2 Kinematic Path Planning Subject to Equal- 
ity Constraints 

In this paper, we consider the problem of finding a continuous path that links the 
specified initial and final configurations while satisfying a set of specified equality and 
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inequality constraints (which may be either holonomic or nonholonomic). To state this 
problem precisely, we first represent a kinematic model as a control system: 

q(t) = u(t) ; g(0) = q 0 (1) 

where q G R n is the configuration variable and u G R n is a pseudo-velocity consid- 
ered as the control input. Note that since this is a path planning problem, the path 
variable, f, is arbitrarily normalized to be within the interval [0,1]. There may also 
be additional holonomic and nonholonomic constraints imposed along the path. If the 
constraints along the path are directly incorporated in the kinematics equation, these 
constraints are treated as hard constraints. If the planner iteratively updates the path 
until the constraints are within a certain tolerance, these constraints are considered as 
soft constraints. We first state the general path planning problem with soft constraints: 

Define u = {u(f) : t G [0,1]} and £ = { q(t ) : t G [0,1]}- Given (1), 
find u G L 2 ([0, l];R n ) such that q(t) satisfies (1) and c(q,u) = 0 where 
c : X 2 ([0, 1]; R n ) X L 2 {[ 0, 1]; R n ) -*Y is a given equality constraint function 
for a specified normed linear space , Y. 

The constraint function c may include physical constraints, such as nonintegrable ve- 
locity constraints on wheels (nonholonomic) and contact constraints of manipulators 
(holonomic), or artificial constraint such as the desired end effector path of a redundant 
manipulator. Pathwise holonomic constraints can be represented as 

c(g,u)(t) = k^qit)) , t G [0, 1] (2) 

and pathwise nonholonomic constraints are of the form 

c(£, u)(t) = k 2 (t, q{t), £(t)) , t G [0, 1] (3) 

and can not be integrated to a constraint only involving the configuration variable. 
In many applications, k 2 is linear in q and invariant in t, i.e., k 2 (t,q(t),q(t)) = 
K 2 (q(t))q(t). A system may also be subject to what is known as second order nonholo- 
nomic constraints (for example, for a manipulator with unactuated joints [47]). 

For the path planning problem with hard constraints, the path constraints in (2) 
and (3) are explicitly removed. In the holonomic case, write (2) in the differential form: 

, dk 2 (t,q) u _ ( 4 ) 

dt 8q 

Denote the Jacobian matrix by J(t,q). If J is a fat matrix (as in the case 

of redundant manipulators) and is nonsingular, then the constraint can be explicitly 
incorporated in the kinematics, resulting in 

j = -/+<», <») 

where J + (t,q) denotes the Moore-Penrose pseudo-inverse of J(t, q), J(t t q) is an arbi- 
trary full rank matrix whose range coincides with the null space of J{t,q ), and v is the 
new effective control variable. 
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In the nonholonomic case, if the constraint is linear in u, then the constraint can 
be eliminated, resulting in a new kinematic equation: 

q = K 2 (q)v (6) 

where v E R m is the new effective control variable (m is the dimension of the null 
space of ^ 2(^)5 assuming it is a constant). 

The general planning problem can now be restated as before except the kinematic 
equation is modified as below to incorporate the path constraints: 

«(<) = &(*. «(*)) + /(*> «(*))«(<) ; g(o) = ?o- (7) 


The drift term, fi(t,q), is zero in the case of nonholonomic constraints linear in q. The 
equality constraint function may involve only the end configuration: 


c(q,u) = q( 1) - q d 


( 8 ) 


where q d is the desired final configuration. In the case of redundant manipulators, joint 
cyclicity implies q d = go- 

For the path planning problem as formulated above with either soft or hard con- 
straints, our approach is based on converting the differential description of the problem 
into an algebraic form. First consider the soft constraint formulation. Recall this case 
involves a very simple kinematic equation given by (1). For a given initial configuration 
qo, q can be related to u via a linear causal map, D, and an initial condition vector, 

So’ ... 

q- Du + q^. (9) 

Define the constraint error as 


2 = ciRu + ^u). 


( 10 ) 


The path planning problem is now reformulated as a nonlinear root finding problem 
for c as a function of u Nonlinear root finding has the reputation being numerically 
challenging [48, §9.6]. However, the situation here is different than the general case 
in that the dimension of the search variable, u, is typically much larger than the 
dimension of the constraint equation. Consequently, there are a very large number 
of roots and finding one of these roots is less difficult than the general root finding 
problem. However, equality constraints alone may not produce physically realizable 
solution; we shall see how this formulation can be extended to inequality constraints 
in the next section. 

To find a u that solves y = 0, our basic strategy is to lift a path in Y that connects 
an initial guess y(0) to zero to a path in £2([0> l]> R m )* Then the end point of the path 
in £ 2 ([0, 1]; R m ) is a solution that satisfies the stated equality constraint c(g, u) - 0. 
To achieve this, we set up the path iteration equation: 


dg dg 

dr dr 


( 11 ) 


where 

G = V g c D_ + Vu£. 
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If G is full rank, or, equivalently, the null space of G* (the adjoint map of G) is zero, 
then any one of the following algorithms can be applied to update u : 

£ = + G( ( 12 ) 

ar ar 

where G + is the Moore-Penrose pseudo-inverse of G, and is any desired convergence 
profile of y in r and j^(0) = y(0), for example, we can choose j^(t) = e- QT y(0), 
for exponential convergence. However, this rate may not correspond to the physical 
convergence rate since (12) may take much longer to solve. In practice, Eq. (12) may 
be solved by an ODE solver or discretized with the time step in r found by line search. 

In the hard constraint formulation and with only the end configuration constraint 
(i.e., c(g,u) = ^(l) - qd) , the equality constraint error can be written as 

y = F(u)-q d (13) 

where F maps the control, u, (for a given < 70 ) to the final configuration. The analytic 
form of F is in general difficult to find, and will not be explicitly required. The system 
(7) is globally controllable if and only if the nonlinear map F is onto, and the system 
is locally controllable around u if and only if ^ uF(u) is a linear onto map. Local 
controllability around u is equivalent to the controllability of the linear time varying 
system obtained after linearizing (7) around q and u. The path planning problem can 
be iteratively solved as before using (12) but with 

G = V„F(tt). 

A sufficient condition for convergence is that G is full rank for all t, or, equivalently, 
the time varying linearized system with u(r ) is controllable. For systems without drift, 
such as nonholonomic systems, it has been shown in [43] that this condition is generic. 

Equation (12) can be solved by an ordinary differential equation (ODE) solver with 
an initial guess of u(0) such that the full rank condition is satisfied. In addition, ^(t) 
can be used to check the accuracy of the solution, y(r), at each r. It can be shown 
that if u(0) is sufficiently close to a solution, then the wellposedness of the IVP (and 
its convergence) is assured. As stated before, the wellposedness of the IVP in the 
hard nonholonomic constraint case has been shown to be generic and there has been 
a surge of recent interest of this technical issue [39]. For the general problem, the 
wellposedness condition remains an important research topic. In all of our simulation 
experience involving equality constraints, this has never been a problem. In the next 
section, when constraints are included through the global exterior penalty function, 
the issue of wellposedness becomes more severe. 

An interesting aspect of this approach is that in (12), £ does not affect the guaran- 
teed convergence rate (specified by a), though it does affect the way ji converges. Since 
the dimension of n is much larger than there is much freedom in £ to affect the even- 
tual convergent solution. For example, £ may be chosen so that additional constraints 
in Ji may be satisfied. We have explored choosing £ via quadratic programming with 
some success, but more research needs to be done in exploiting this degree of freedom. 

There are other possible choices of 34 , f° r example, in the soft constraint case. 

§ = G‘^ + G* ( 14 ) 

OX OT 
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du 

dr 


(15) 


G'i 




dw d 

dr 


2 ^f + <?e 


^ = i|NI 2 


The trade-offs between computational load and convergence speeds between these 
schemes are yet to be fully explored. 

The initial guess tt(0) will clearly affect the convergence. We have not extensively 
explored an intelligent procedure (perhaps based on past experience) for this selection. 
In all of our examples, the initial conditions are simply chosen to avoid the rank 
deficiency of G. 

In the algorithm described above, it is critical to be able to calculate G efficiently. 
Since G relates an infinitesimal variation Su to a corresponding infinitesimal variation 
Sc, it can be computed based on the linearized kinematic equation (and the constraint 
equation in the soft constraint case). This procedure is described in detail in [49], and 
has been used in all of our examples. 


3 Kinematic Path Planning Subject to Inequal- 
ity Constraints 

In the previous section, only equality constraints are considered. However, for realistic 
problems, there are also many inequality constraints which need to be enforced. For 
problems involving equality constraints only, the obtained path is frequently undesir- 
able. For example, in the case of a tractor with twin-trailers, the front wheel angles of 
the tractor may go through an unrealistically large range motion, the trailers can be- 
come jack-knifed, and the entire vehicle may go through several complete revolutions; 
in the redundant arm case, the joints can violate joint limits, and there may be self 
collision or collision with other objects in the workspace. Clearly, if this algorithm is 
to be used in practical applications, additional constraints must be incorporated. In 
this section, we present an approach similar to the exterior penalty function method of 
[45] and which has so far been shown to be very effective in addressing the inequality 
constraint issue. 

There are three basic approaches to address the inequality constraints: 

• Convert the inequality constraints into equality constraints by defining a function 
that is zero when the constraint is satisfied and non-zero when the constraint is not 
satisfied. Once this is done, the zero finding approach of the previous section can be 
applied to find a zero solution which corresponds to a feasible path. 

• If the feasible region is convex polyhedral, then linear programming or quadratic 
programming can be used to select the free variable £ in (12) without affecting the 
convergence of the equality constraints. 

• Pre-warp the vector with a multiplicative potential function <j>{q) which is zero on 
the boundary of the inequality constraint and positive inside the feasible region. It has 
recently been pointed out in [44] that as long as the initial path q( 0) lies strictly within 
the feasible region, then the same algorithm (12) can be applied to iteratively move y 
to zero. 

We have extensive experience with the first two approaches. The first approach has 
proven to be very effective,, though there are currently no (genericity) results on the 
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full rank gradient condition as in the equality constraint case. The second approach is 
computationally intensive due to the need to repeatedly solve a constrained optimiza- 
tion problem. It is perhaps best used as a local optimization refinement in conjunction 
with the first approach. The third approach has just recently been suggested and is 
computationally untried; we intend to pursue this direction in the proposed research. 
For the remainder of this section, we will elaborate on the exterior penalty function 
approach. 

Suppose the feasible region is defined by a set of p inequalities: 

5i(2,«) <0 (16) 

where g\ is nonlinear and assumed to be as smooth as needed, and < is interpreted 
as a component-wise relationship. These constraints may be directly placed on the 
configuration variables such as joint limits or on other variables such as the end effector 
of a manipulator, vehicle boundaries, etc. Without loss of generality, we shall consider 
the hard constraint case for the discussion below. Let the relationship between ^ and u 
be denoted by q = F(u). By substituting this relationship, the constraint inequalities 
are transformed to 

g(F(u),u) <0 (17) 

or more simply expressed as 

< 0. (18) 

Suppose now we wish to constrain the states of the system to stay within the feasible 
region defined by the task space. We then define a penalty function corresponding to 
the i th system state as: 


_ V / 7 j« 0 j,U) 9 ji(u) > 0 

' £il 0 sM < o 


(19) 


where A > 1, gj^s correspond to the constraint applied to state i at discrete time j and 
7 i» s are constant weights. This function is nonzero when state i violates the constraints 
at any time along the path, and is identically zero only when all constraints are satisfied. 
For this penalty function, the penalty imposed depends upon the constraint violation 
raised to the A power. Therefore, the penalty function approaches infinity as the system 
makes ever increasing incursion out of the feasible region. Of course many other choices 
for the constraint function are also possible. For example, another penalty function 
which we use frequently is the following: 


z, = rf *.(*)> o 

h\ o #*(■) < o 


( 20 ) 


where r, > 1 and A > 0. This penalty function is bounded by the sum of the 7 j,’s, and 
is therefore in some sense less harsh than the previous penalty function. In practice 
we have seen that the second penalty function gives faster convergence than the first, 
in cases where the vehicle makes large excursions outside of the feasible region. This 
advantage in convergence stems from the fact that is smaller for the second 

constraint function than for the first when large constraint violations occur, causing 


12 


the IVP to be more wellposed for the second constraint function. For either penalty 
function, the composite constraint vector is defined as z = [ z\ ... z n T . 

We now have a situation similar to the unconstrained case: Find u so that y = 0 
and z- 0. The same approaches described in Section 2 can now be applied. Using the 
hard constraint formulation, the new differential equation in u now becomes: 


dy^r) 

dr 


Gt(u(r)) 


h 

*d 


+ Gi(u(t))((t) 


( 21 ) 


where, 



( 22 ) 


and m and z d describe the desired path in r that links the initial error to zero. The 
gradient of z( u) can be easily obtained via chain rule and the linearized kinematic 
equation [49]. The convergence of this algorithm now requires the full rank condition 
on G\ in each iteration. Note that in (21), there is again a free parameter £(t). The 
extra degrees of freedom offered by f(r) may be used to satisfy additional optimality 
considerations. 

The penalty function formulation we have presented allows for a wide range of 
constraint types. Up to this point we have considered constraints which apply directly 
to the configuration variable of the system, but the formulation allows constraints to 
be applied to non-configuration variables as well. In any practical application it is 
not enough just to constrain the origin of the system inside a feasible region. The 
boundaries of the system must stay inside as well. For instance, for a single planar 
body, we relate any point p on the boundary to the configuration variable q by the 
nonlinear transformation: 


Px 

Py 


?1 

+ 

cos(g 4 ) -sin(g 4 ) 

r x 

q 2 

_ sin(g 4 ) cos (q 4 ) 

. r y . 


(23) 


where r x and r y are the x and y positions of the boundary point in the body frame. 
The same penalty function which applied to q before now applies to p, gj{p(qj («)))• 
A similar transformation can be obtained to relate boundary points of multi-bodied 
systems to the configuration variable. 

The penalty function formulation also allows for a wide range of feasible region 
types. For instance, when the feasible region is polyhedral, g takes on particularly 
simple form: 

g(q) = Ag-b. 

We have looked at this type of constraint extensively. In fact, this constraint type is 
always used when limiting wheel or jackknife angles on wheeled vehicles. This type of 
constraint is simple to use but the feasible region is necessarily convex: 

The penalty function can be used to enforce both convex and non-con vex, non- 
polyhedral constraints. For instance, supposed it is desired to drive the system while 
keeping boundary points outside of an circular region. The constraint g at each time 
j can then be expressed as: 


9j(p ) = ~(Px - X 0 ) 2 - ( Py - yo) 2 + r 2 
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where xo and yo are constants, r is the radius, and p is a boundary point as described in 
(23). We have had much success in using this formulation with cars and tractor-trailer 
vehicles and present an example in the next section. 

In the preceding discussion, system constraints are expressed as analytic functions 
of the configuration variable. 'But again, in practical situations it may not be possible 
to adequately represent the constraints by analytic functions. In particular, for highly 
unstructured environments it may require great effort to find an appropriate analytic 
function, and once found, it may be computationally intensive to apply it. Therefore, 
we propose a method based upon the contour map built up from the task space obstacle 
boundaries. First, suppose that the feasible region is defined by a set of p inequalities 
of the form: 

92(Z(u)) < 0 (24) 

where g 2 is nonlinear, smooth as needed, and < is component-wise. The gradient of 
02 with respect to u is then obtained by applying the chain rule, 

Vu52 = [V lt Z] T V i 0 2 . (25) 

Note that in this equation, is independent of the constraints and that 

depends solely upon the constraints. Therefore, rather than compute g 2 and V q g 2 
explicitly using analytic functions, it is possible to use a lookup table to compute these 
values. For < 72 , the lookup table contains the cost for each point in a grid covering the 
task space. For the table contains the gradient components for each point on the 
grid. As before, this formulation will also work for non-configuration variables. Since 
the exterior penalty function approach is used, the contour map is zero everywhere 
inside the feasible region and non— zero outside. We have just begun using this method 
for path planning and have so far had good success. An example of the contour map 
method applied to a car is given in the next section. 

In this, and the previous sections, we have used the notation u to denote the set 
of control inputs applied to the system at each point in time along the path. For 
the continuous case this set contains an infinite number of elements. Therefore u in 
vector form would be an infinite dimensional vector. However, in order to implement 
the algorithm on a digital computer u needs to have finite length. The first obvious 
step would be to discretize u(t) using the standard basis. We have used this approach 
extensively and have obtained good results from it. Another approach we have tried is 
to use the first N elements of the Fourier basis to approximate u(t). In this formulation, 
for each discrete time j we represent the control input as: 

uj - ( 26 ) 

where $ is a matrix containing the Fourier basis elements and A is the constant vector 
of Fourier coefficients. In this formulation, A uniquely describes the control over all 
time, since it is independent of time. Now for a given initial configuration g can be 
related to A by a causal map, 

* = Zi(A>. (27) 

Furthermore, substituting for £ and A for U. in all previous equations, we see 
that all previous methods still work. The full advantages and disadvantages in using 


the Fourier basis formulation rather than the discrete basis is unknown at present 
However, in the case of driving a car around a circle, we have observed that the discrete 
basis formulation has great difficulty arriving at a solution whereas the Fourier basis 
formulation solves the problem with relative ease. 


4 Examples 


We have applied the algorithm to a large number of computer simulation examples 
and experimentally to a mobile robot. For illustration purposes, we include several 
examples involving three different wheeled vehicles and a point-to-point path planning 
of a 4R planar arm and a 9DOF arm like the one in our lab (see next section on a 
fuller description of this arm). 

In the first example, we consider the parallel parking of a double tractor-trailer. A 
true double tractor-trailer, like those seen driving on the highways, actually consists 
of a tractor with three trailers: two long trailers connected by a comparatively short 
trailer. The short trailer, or dolly, makes any backing-up situation extremely difficult 
in that small backward motions can produce large jackknife angles between the dolly 
and the trailers. In fact, it is known that with a human operator one must, in general, 
disconnect the end trailer and dolly before any backward motions are attempted. 

In the example below, shown in Fig. 1-2 below, the front wheels are constrained 
to ±35° and the jackknife angles (the jackknife angle is defined as the angle between 
the center line of one trailer and the center line of an adjacent trailer) are constrained 
to ±50°. The limits on the jackknife angles are sufficient to ensure that none of 
the trailers collide. In this example the algorithm required about four hundred fifty 
evaluations of the right hand side of the differential equation to achieve a tolerance on 
the potential fields of 0.01. The penalty function method used here is simply the first 
convex, polyhedral method mentioned in the previous section. 


S 


Figure 1: Parallel Parking Path of a Double Tractor-Trailer with Constraints 

The next example is presented to demonstrate the versatility of the algorithm with 
regard to the type of constraint applied. In this example, it is desired to drive a 
tractor-trailer vehicle around a circular obstacle while remaining within a rectangular 
region limiting motion along the z-axis. To further complicate the example, two other 
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Figure 2: Parallel Parking Path of a Double Tractor-Trailer with Constraints 


circular obstacles are included above and below the first one. In addition, these task 
space constraints are applied to fourteen different non-configuration variable points 
around both the boundary of the tractor and the boundary of the trailer. The steering 
angle is constrained to ±20°, and the jackknife angle is constrained to ±50°. We 
therefore have in this one example, all of the different types of constraints mentioned 
in the previous section: convex polyhedral (steering, jackknife, and x-position), non- 
convex non-polyhedral (circular obstacles), and non-configuration variable (boundary 
point) constraints. The Fourier basis representation of the control input, with thirty 
three basis elements, was also used in this example. Using a steepest descent method 
with a Golden section search, the algorithm took fifty four evaluations of the right hand 
side to converge. The initial path in this case, drives straight through the obstacle. 
The final result is as shown below in Fig. 3. 


solkJ-xl. da»h-x2. dot*x3, djsh-dcc*x4, + nj-x4 
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Figure 3: Tractor-Trailer Example with Non-Convex Task Space Constraints 
The next simulation example uses the contour map method to calculate the penalty 
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function. In this example, it is desired to drive a car around a circular obstacle while 
remaining inside of a circular region. Seven points around the boundary of the car are 
used^for non-configuration variable constraints. The front wheels are constrained to 
±15°. The Fourier basis was once again used to represent the control input along the 
path. This example took about one hundred evaluations of the right hand side using 
a discretized version of (21) with Golden section search to find the optimal step size 
The initial path drives straight through the obstacle. Below, the first set of plots, in 
Fig. 4, show the car and path in the context of the task space. The second set of plots, 
in Fig. 5 show the contour map used in this example. 
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Figure 4: Path Planning with Contour Map Method 



Figure 5: Plots of Exterior Penalty Function and Contour Map 
In the final wheeled vehicle example, we present experimental results of applying the 
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path planning algorithm to an actual mobile robot. The robot is actually a one quarter 
scale model car, complete with a passive suspension system, on-board IBM compatible 
486 33-MHz computer, and wheel encoders for use in dead-reckoning feedback. The 
desired path was a simple parallel parking path with constraints only on the steering 
angle of ±15° which is the physical limit of the car. The purpose of this experiment 
was to determine if a real robot can actually follow a path generated by the planning 
algorithm. In this example, the path was generated, and then a velocity profile was 
imposed. A simple, Lyapunov function based, feedback controller was used to try and 
track the path. Although there is some error in tracking the path, it is obvious from 
the plots below, in Fig. 6-7 that the path is such that the real car can follow it. The 
tracking errors are due in part to error in the dead-reckoning feedback and to the 
controller used. 


solkJ*xtual. dotted-desired, units in inches 



-20 0 20 40 60 80 100 

s-poution 


Figure 6: Parallel Parking of Catmobile: Experimental vs. Simulation Path 

Various examples of redundant manipulators are presented which incorporate both 
joint and task space inequality constraints. Key features of this approach include vari- 
ous possible goal task specifications (from end-point only to entire path specification), 
joint path cyclicity, and robustness to manipulator singularities. Inequality constraints 
are handled by the global penalty functions and a linear inequality set description as 
shown in Sections 3. Simulations were performed using Matlab with computationally 
intensive routines coded as Matlab-callable C routines (cmex files). Note that the end 
effector constraint is treated as a soft constraint here as compared to the nonholonomic 
constraint considered above which is treated as a hard constraint. 

The first example illustrates the joint sequence for a 3R planar arm required to trace 
out the tip path shown while remaining within a set of task space boundaries (Fig. 8). 
The path shown in Fig. 8 includes an intermediate joint vector in which links 2 and 3 
become aligned - equivalent to a pose switch for spatial arms. Local planning methods 
typically encounter difficulty in handling pose changes since they correspond to arm 
singularities, and the arm Jacobian losing rank. The present algorithm executes the 
pose switch smoothly, while tracking the desired tip path. In the next situation (Fig. 9) 
a 4R arm is required to traverse to specified coordinates Xj . The obstacles are chosen 
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Figure 7: Parallel Parking of Catmobile: Experimental vs. Simulation State Trajectories 

to provide a relatively narrow opening through which the arm must pass. Typically, 
this type of scenario is very challenging for planners based on purely local potential 
field formulations due to a local minimum formed in the space between Obstacles 1 
and 2. By iterating from an intermediate path sequence generated for the goal end- 
point without consideration of the obstacles, it is possible to warp the intermediate 
path to meet all the constraints. It should also be noted that the final path cannot be 
accomplished without switching the pose along the way. 



Figure 8: Cartesian Tip Path Figure 9: Obstacle Avoidance Problem 


The ability to incorporate joint path cyclicity as an equality constraint is illustrated 
by an example where a 3R arm must trace out the boundaries of a Cartesian square 
path, within a restricted workspace (Fig. 10). Path planning for cooperating arms ma- 
nipulating a common payload imposes an additional kinematic closure constraint on 
the arm tips. Fig. 11 shows the path sequence generated for a pair of 3R planar manip- 
ulators required to move the connecting linkage from its initial position of ( — 1, 2, 0°) to 
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a goal position of (2,2,45°). In the planar examples shown, the apparent collisions of 
the links with the obstacle boundaries arise out of the fact that checking is performed 
for the tips of each link only. A more complete collision detection procedure is being 
developed for spatial arm planning. 



Figure 10. 3R Cyclic Path Figure 11: Cooperating Arm Problem 

The remaining examples illustrate applications to a spatial redundant 9 DOF arm. 
This consists of a PUMA 600 manipulator mounted on a platform providing additional 
Unear, rotate and tilt joints. (See [50] for details of the CIRSSE dual-arm robotic 
testbed.) Fig. 12 shows the output sequence for joint path end-point planning to a 
specified task space position/orientation in the presence of an obstacle. Fig. 13 shows 
the path sequence generated for a path foUowing task incorporating a straight Une 
translation coupled with a rotation of 180° about the tip Y axis. Because of the 
manipulator joint Umits, this can only be accompUshed by switching the PUMA’s pose 
from elbow-up to elbow-down at some intermediate point. As in the planar examples, the 
present algorithm accompUshes a smooth transition between these arm configurations. 
The final example (Fig. 14) illustrates a cycUc joint path sequence computed for the 
arm to follow a circular task space path while maintaining a fixed tool orientation 
(perpendicular to the plane of the circle). 

For the cases shown, the planar examples required from seven to ten iterations to 
reduce the task error and meet all the constraints, while the spatial examples ranged 
from ten to twenty-five iterations. Discretization levels ranged from N = 10 to 40 
in size. For the 9 DOF redundant arm planning scenario, this represents a small 
computational load for the Sun SPARC-station used, workstations. For problems of 
larger size, this technique can be applied recursively to the desired discretization level, 
thereby keeping the individual iteration array operations small. The ability to apply 
the algorithm to the global problem (even at the coarsest resolution) is essential for 
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Figure 12: 9 DOF Path Planning Figure 13: 9 DOF Path Following 

developing globally feasible solutions to the constrained path planning problem. 


5 Conclusions 

Kinematic path planning for robots is a key step in their effective utilization on earth 
and in space. The various types of constraints such as holonomic and nonholonomic, 
and equality and inequality constraints pose particular challenges. In this paper, we 
consider a novel and promising method which warps the entire path until all constraints 
are satisfied. The main approach is to convert the differential (local) kinematic rela- 
tionship to an algebraic (global) equation. With emphasis placed on feasibility rather 
than optimality, we obtain an initial value problem rather than two-point boundary 
value problem typically arisen in a global optimization approach. This formulation is 
general enough to include both redundant manipulators and nonholonomic systems, 
and combination of the two; these topics are traditionally treated separately due to 
their unique properties when local algorithms are applied. The inequality constraints 
are handled through a global exterior penalty function method, which allows for non- 
polyhedral constraints as well as constraints on non-configuration variables. For the 
future research, we will address the following fundamental issues related to this promis- 
ing algorithm: 

1. Develop conditions of convergence of the algorithm based on the wellposedness of 
the initial value problem. Also develop strategy to proceed when the algorithm 
fails to converge due to the rank deficiency of G in (12). 

2. Develop an algorithm to adaptively adjust the path discretization step size based 
on the local discretization error. 

3. Incorporate optimality as a secondary criterion. 

4. Improve the sensitivity and robustness of the algorithm with respect to imperfec- 
tion in the kinematic model. 

5. Implement the planner on various experimental platforms. 
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Figure 14: Spatial 9 DOF Cyclic Cartesian Path 

Acknowledgment 

This work is supported in part by the Center for Intelligent Robotic Systems for Space 
Exploration and New York State Center for Advanced Technology for Automation and 
Robotics at Rensselaer Polytechnic Institute. 


References 

[1] Lozano-Perez T. and Taylor R.H. Geometric issues in planning robot tasks. In 
Problems of Robotics. MIT Press, Cambridge, Mass., 1989. 

[2] J.C. Latombe. Robot Motion Planning. Kluwer International Series in Engineer- 
ing and Computer Science: Robotics: Vision, Manipulation and Sensors. Kluwer 
Academic Publishers, 1991. 

[3] Hwang Y.K. and Ahuja N. Gross motion planning - a survey. ACM Computing 
Surveys, 24(3):219-291, 1992. 

[4] Lozano-Perez T. and Wesley M.A. An algorithm for planning collision free paths 
among polyhedral obstacles. Communications of the ACM , 22(10):560-570, Oc- 
tober 1979. 

[5] Mazer E. Lozano-Perez T., Jones J.L and O’Donnell P.A. Task-level planning of 
pick-and-place robot motions. Computer, 22(3):21-29, March 1989. 


22 


[6] Weaver J.M. and Derby S.J. A divide-and-conquer method for planning collision- 
free paths for cooperating robots. In AIAA Space Programs and Technologies 
Conference, Paper AIAA-92-1722, Huntsville, AL, 1992. 

[7] Whitney D.E. Resolved motion rate control of manipulators and human prosthe- 
ses. IEEE Transactions on Man-Machine Systems, MMS-10(2):47-53, 1969. 

[8] Khatib 0. Dynamic control of manipulators in operational space. In 6th IFToMM 
Congress on Machines and Mechanisms, New Delhi, 1983. 

[9] Vukobratovic M. and Kircanski M. A dynamic approach to nominal trajectory 
synthesis for redundant manipulators. IEEE Transactions on Systems, Man and 
Cybernetics, SMC- 14, 1984. 

[10] Yoshikawa T. Manipulability of robot mechanisms. International Journal of 
Robotics Research, 4(2):3-9, March 1985. 

[11] Hollerbach J. Baillieul J. and Brockett R. Programming and control of kinemat- 
ically redundant manipulators. In Proc. 24th IEEE Conf. Decision and Control, 
pages 768-774, Las Vegas, NV, December 1984. 

[12] Hollerbach J.M. and Suh K.C. Redundancy resolution of manipulators through 
torque optimization. In Proc. 1985 IEEE Robotics and Automation Conference, 
pages 308-315, St. Louis, MO, March 1985. 

[13] B ailli eul J. Martin D.P. and Hollerbach J.M. Resolution of kinematic redun- 
dancy using optimization techniques. IEEE Transactions on Automatic Control, 
5(4):529-533, August 1989. 

[14] A. De Luca. Redundant Robots: Introduction to Chapter IX. In M. Spong and 
F. Lewis, editors, Robotics. IEEE Press, 1992. 

[15] D.P. Martin, J. Baillieul, and J.M. Hollerbach. Resolution of kinematic redun- 
dancy using optimization techniques. IEEE Transactions on Robotics and Au- 
tomation, 5(4), August 1989. 

[16] Y. Nakamura. Advanced Robotics: Redundancy and Optimization. Addison— 
Wesley Series in Electrical and Computer Engineering: Control Engineering. 
Addison-Wesley Publishing Company, 1991. 

[17] Suh K.C. and Hollerbach J.M. Local versus global torque optimization of redun- 
dant manipulators. In Proc. 1987 IEEE Robotics and Automation Conference, 
pages 619-624, Raleigh, NC, March 1987. 

[18] Kazerounian K. and Nedungadi A. Redundancy resolution of robotic manipulators 
at the acceleration level. In The 7th World Congress on the Theory of Machines 
and Mechanisms, pages 1207-1211, Sevilla, Spain, September 1987. 

[19] Wang Z. and Kazerounian K. A general formulation for redundancy resolution 
of serial manipulators. In Proc. 1990 ASME Design Technical Conference, pages 
373-379, Proc. 1990 ASME Design Technical Conference, September 1990. 

[20] P. Jacobs and J. Canny. Robust motion planning for mobile robots. In 1991 IEEE 
R&A Workshop on Nonholonomic Motion Planning, Sacramento, CA, April 1991. 

[21] J. Barraquand, B. Langlois, and J.-C. Latombe. Numerical potential held tech- 
niques for robot path planning. In Proc. IEEE 5th Int. Conf. on Advanced 
Robotics, pages 1012-1027, Pisa, Italy, June 1991. 


23 


[22] P. Jacob, J.-P. Laumond, and M. Taix. A complete iterative motion planner for 
a car-like robot. Joumees Geometrie Algorithmique, 1990. 

[23] L. Dorst, I. Mandhyan, and K. Trovato. The geometrical representation of path 
planning problems. Technical report, Philips Laboratories North American Philips 
Corp., Briarcliff Manor, New York, 1991. 

[24] B. Mirtich and J. Canny. Using skeletons for nonholonomic path planning among 
obstacles. In Proc. 1992 IEEE Robotics and Automation Conference, pages 2533- 
2540, Nice, France, May 1992. 

[25] T. Hague and S. Cameron. Motion planning for nonholonomic industrial robot 
vehicles. In IEEE/RSJ International Workshop on Intelligent Robots and Systems, 
pages 1275-1280, Osaka, Japan, November 1991. 

[26] T. Fraichard and C. Laugier. On line reactive planning for nonholonomic mobile 
in a dynamic world. In Proc. 1991 IEEE Robotics and Automation Conference, 
pages 432-437, Sacramento, CA, April 1991. 

[27] R.W. Brockett. Control theory and singular Riemannian geometry. In New Di- 
rections in Applied Mathematics, pages 11-27. Springer- Verlag, New York, 1981. 

[28] Z.X. Li and J. Canny. Motion of two rigid bodies with rolling constraint. IEEE 
Transactions on Robotics and Automation, 6(l):62-72, February 1990. 

[29] G. Lafferriere and H.J. Sussmann. Motion planning for controllable systems with- 
out drift: A preliminary report. Technical Report SYCON-90-04, Rutgers Center 
for Systems and Control, June 1990. 

[30] R.M. Murray and S.S. Sastry. Steering nonholonomic systems using sinusoids. In 
Proc. 29th IEEE Conference on Decision and Control, pages 2097-2101, Honolulu, 
HI, 1990. 

[31] Z. Li and L. Gurvits. Theory and applications of nonholonomic motion plan- 
ning. Technical report, New York University, Courant Institute of Mathematical 
Sciences, July 1990. 

[32] Y. Nakamura and R. Mukherjee. Nonholonomic motion planning of space robots. 
In Proc. 1990 IEEE Robotics and Automation Conference, Cincinnati, OH, 1990. 

[33] Z. Vafa and S. Dubowsky. On the dynamics of manipulators in space using the 
virtual manipulator approach. In Proc. 1987 IEEE Robotics and Automation 
Conference, Raleigh, NC, March 1987. 

[34] A.M. Bloch, N.H. McClamroch, and M. Reyhanoglu. Controllability and sta- 
bilizability properties of a nonholonomic control system. In Proc. 29th IEEE 
Conference on Decision and Control, pages 1312-1314, Honolulu, HI, 1990. 

[35] B. d’Andrea Novel, G. Bastin, and G. Campion. Modelling and control of non- 
holonomic wheeled mobile robots. In Proc. 1991 IEEE Robotics and Automation 
Conference, pages 1130-1135, Proc. 1991 IEEE Robotics and Automation Con- 
ference, April 1991. 

[36] Y. Kanayama, Y. Kimura, F. Miyazaki, and T. Noguchi. A stable tracking con- 
trol method for an autonomous mobile robot. In Proc. 1990 IEEE Robotics and 
Automation Conference, pages 384-389, Cincinnati, OH, May 1990. 


2b 


• * 


[37] C. Samson and K. Ait-Abderrahim. Feedback control of a nonholonomic wheeled 
cart in cartesian space. In Proc. 1991 IEEE Robotics and Automation Conference , 
Sacramento, CA, April 1991. 

[38] C. Samson and K. Ait-Abderrahim. Feedback stabilization of a nonholonomic 
wheeled mobile robot. In IEEE/RSJ International Workshop on Intelligent Robots 
and Systems, pages 1242-1247, Osaka, Japan, November 1991. 

[39] H.J. Sussmann and Y. Chitour. A continuation method for nonholonomic path 
finding problem. IMA Workshop on Robotics, January 1993. 

[40] A. Divelbiss and J.T. Wen. A perturbation refinement method for nonholonomic 
motion planning. In Proc. 1992 American Control Conference, Chicago, IL, June 
1992. 

[41] A. Divelbiss and J.T. Wen. Nonholonomic motion planning with constraint han- 
dling: Application to multiple-trailer vehicles. In Proc. 31th IEEE Conference on 
Decision and Control, Tucson, AZ, December 1992. 

[42] S. Seereeram and J.T. Wen. A global approach to path planning for redundant 
manipulators. Accepted for publication in IEEE Trans, on Robotics and Automa- 
tion, 1993. 

[43] E.D. Sontag and Y. Lin. Gradient techniques for systems with no drift. In Pro- 
ceedings of Conference in Signals and Systems, 1992. 

[44] E.D. Sontag. Non-singular trajectories, path planning, and time-varying feedback 
for analytic systems without drift. IMA Workshop on Robotics, January 1993. 

[45] Dorny C.N. A Vector Space Approach to Models and Optimization. R. E. Krieger 
Publishing Co., Florida, 1986. 

[46] D.G. Luenberger. Linear and Nonlinear Programming. Addison- Wesley, 2 edition, 
1984. 

[47] G. Oriolo and Y. Nakamura. Free-joint manipulators: Motion control under 
second-order nonholonomic constraints. In IEEE/RSJ International Workshop 
on Intelligent Robots and Systems, Osaka, Japan, November 1991. 

[48] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. Numerical 
Recipes: The Art of Scientific Computing. Cambridge University Press, Cam- 
bridge, U.K., 1986. 

[49] A. Divelbiss and J.T. Wen. A perturbation refinement method for nonholonomic 
motion planning. CAT report #10, Rensselaer Polytechnic Institute, March 1992. 

[50] J.F. Ill Watson, D.R. Lefebvre, S.H. Murphy A.A. Desrochers, and K.R. Field- 
house. Testbed for cooperative robotic manipulators. In A.A. Desrochers, editor, 
Intelligent Robotic Systems for Space Exploration. Kluwer Academic Publishers, 
1992. 


25 


Application of the CIRSSE Cooperating 
Robot Path Planner to the NASA 
Langley Truss Assembly Problem 

Jonathan M. Weaver and Stephen J. Derby 
Rensselaer Polytechnic Institute 
Troy, New York 



PJWdOiNG PAGE BLANK NOT FILMED 


N94- 26280 


/ c ) c (^b <P/7 V 


Application of the CIRSSE Cooperating Robot Path Planner 
to the NASA Langley Truss Assembly Problem 


Jonathan M. Weaver Stephen J. Derby 

Research Assistant Associate Professor 

Center for Intelligent Robotic Systems for Space Exploration 
Department of Mechanical Engineering, Aeronautical Engineering, and Mechanics 
Rensselaer Polytechnic Institute, Troy, New York 12180-3590 


Abstract 

A method for autonomously planning collision /rec 
paths for two cooperating robots in a static environ- 
ment has been developed at CIRSSE. The method uti- 
lizes a divide- and- conquer type of heuristic and in- 
volves non- exhaustive mapping of configuration space. 
While there is no guarantee of finding a solution , the 
planner has been successfully applied to a variety of 
problems including two cooperating 9 dof robots. 

Although developed primarily for cooperating robots , 
the method is also applicable to single robot path plan- 
ning problems. A single 6 dof version of the planner 
has been implemented for the truss assembly task at 
NASA Langley f s Automated Structural Assembly Lab 
(ASAL). The results indicate that the planner could 
be very useful in addressing the ASAL path planning 
problem and that further work along these lines is war- 
ranted. 

1 Introduction 

The robot path planning problem involves deter- 
mining if a continuous and obstacle avoiding path ex- 
ists between start and goal positions, and, if so, to 
find such a path. The complexity of the path plan- 
ning problem has been shown to be exponential in 
the number of dof [1, 2]. A review of the many path 
planning techniques is beyond the scope of this paper. 
Reference [3] presents a recent survey paper on the 
subject. 

A method has been developed at Rensselaer’s Cen- 
ter for Intelligent Robotic Systems for Space Ex- 
ploration (CIRSSE) to autonomously plan collision 
free paths for two robots working cooperatively in 
a known, static environment [4, 5, 6). Cooperation 
refers to the case whereby both robots simultaneously 
grasp and manipulate a common, rigid, payload. The 
planner is based around a divide-and-conquer heuristic 
aimed at traversing c-space while performing selective 
mapping on an as-needed basis. This path planner 
has been applied to the CIRSSE testbed. The testbed 
consists of two 9 dof robots, each of which consists of a 
3 dof platform and a 6 dof Puma. A sample path found 
by the cooperating 9 dof planner is shown in Figure 1. 
This example required approximately 10 minutes so- 
lution time on a SparcStation 1. 


Although developed primarily for the cooperat- 
ing robot case, the c-space traversal heuristic around 
which the planner is based may also be applied 
to single robot path planning problems. This pa- 
per discusses a single arm version of the planner 
which was implemented for the truss assembly task 
at NASA Langley’s Automated Structural Assembly 
Lab (ASAL). The purpose of the implementation was 
to assess the potential usefulness of the planner for the 
ASAL path planning problem. 

The ASAL path planning problem is described in 
Section 2. Our path planning strategy is discussed in 
Section 3. Sections 4 and 5 present some implemen- 
tation details and results for application of the plan- 
ner to the ASAL path planning problem. Section 6 
presents some conclusions and areas for future work. 

2 Problem Statement 

A CimStation model of NASA Langley's ASAL is 
shown in Figure 2 (CimStation model provided by 
NASA Langley). The system consists of a 6 dof Merlin 
robot, shown in Figure 3, mounted to a xy-positioning 
table (referred to as the carriage), and a turntable. 
The turntable includes a triangular platform which 
can rotate around a vertical axis through its center. 
The Merlin robot is kinematically similar to a Puma. 
The objective of the ASAL is to assemble truss struc- 
tures consisting of 102 2 meter long struts. Such a 
truss is illustrated in Figure 4. The truss is assem- 
bled upon the turntable of the ASAL by positioning 
the carriage and the turntable such that the Merlin 
may take each strut from a canister near the base of 
the Merlin and install it in its final position in the 
assembly. 

The ASAL path planning problem as addressed 
herein is defined as follows: Given a carriage and 
turntable position for each strut, determine a suitable 
path for the Merlin to safely move the robot and its 
payload from a start position to a prescribed goal po- 
sition. The start position is above the canister holding 
the as-vet unassembled struts. The goal position for 
each strut is taken as 10 cm from the final position in 
the negative of the approach direction. The assembly 
sequence is as specified by NASA Langley. 

It is assumed that feasible and collision free start 
and goal joint configurations of the robot are known. 
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(a) Start Position 


(b) 



Figure 1: Sample Results for Cooperating 9 DOF Robots 
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Figure 4: 102 Strut. Truss Structure 
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3 Strategy 

Like the single robot planner presented by 
Dupont [7], the principle strategy of our planner is to 
minimize the computationally expensive mapping of 
configuration space by performing mapping on an as 
required basis. The approach is based around a divide- 
and-conquer style heuristic for traversing through 
c-space. Computationally expensive precomputations 
and exhaustive c-space mappings are avoided. The 
approach is applicable regardless of the number and 
type of joints in the robot and for any number of ob- 
stacles in the workspace. A siring tightening algorithm 
may be applied to modify any safe path found by the 
planner into a more efficient one, where efficiency is 
measured by joint space trajectory length. 

The path planning method involves first attempt- 
ing to traverse a c-space vector from the start to the 
goal of one of the robots. If this vector passes through 
unsafe space, the hyperspace orthogonal to and bi- 
secting the unsafe segment of the vector is systemat- 
ically searched to identify an intermediate goal point 
for consideration as a via point. An attempt is made 
to traverse from the last safe point to the intermediate 
goal point. This process is repeated as necessary un- 
til the attempted traversal to the newest intermediate 
goal point is entirely safe. At that point, progres- 
sion is attempted toward all previous guide points in 
the opposite order in which they were found, where 
guide points include not only previous intermediate 
goal points but also the safe points found on the goal 
end of each unsafe region which invoked a search. 
When progression to a particular guide point is not 
entirely safe, that point is permanently dismissed and 
progression is attempted toward the next guide point 
in the specified sequence. The progression continues 
until an attempt has been made to progress to the 
global goal point. If that attempted progression is not 
entirely successful the overall process is repeated until 
the global goal point has been safely traversed to. 

In 2D, the hyperspace orthogonal to an unsafe vec- 
tor (the space which the heuristic searches) is simply 
a line. For 2D problems, the initial search is per- 
formed equally in both directions until a safe point 
is found. Subsequent searches will first exhaustively 
search in the direction which has a component in the 
previously successful search direction. Only when no 
safe point can be found in that direction will the other 
direction be considered. A 2D example of the c-space 
traversal heuristic is shown in Figure 5. This example 
involves non-disjoint safe space and requires multiple 
searches. More 2D examples and a vector description 
of the heuristic may be found in [6]. 

In the general nD case, the search space will be 
n ~ 1 dimensional. In this case, several approaches 
were considered for computing search directions. The 
most effective method found involves considering all 
combinations of ±1 and 0 (except all zeros) times a 
set of orthogonal basis vectors for the subspace. This 

yields 3 n ^ — 1 search directions for an n dof prob- 
lem. The following vectors may be calculated in the 
sequence shown and then normalized to vield one such 
orthogonal basis: 


Bi = (1^1,0 ,0) 

b 2 = ( 6 l 1 »P2’ /l 2»° -* 0 ) 

B 3 = ( & 2 1 « 6 2 2 ’P3« a 3> 0 -- m0) (1) 


B n— 1 = (^n— 2i>^n — 22»* ’ 2 n _2’ 

Pn — 1* ^n-l) 

where the p,* are chosen so that the B,* and B,_ i are 
orthogonal, then the h ^ are chosen so that the B, lie 
in the search hyperplane. 

Initial searches favor all directions equally, whereas 
subsequent searches sort the i directions S 2 - into g 
equal breadth bins by the following rule: 


Sj € bin(j) if ^ < 


d Pi ~ d Pmin 


dpmax — dp, 


min 


- , < 2 > 


where dpj is the dot product of S z * with the previously 
successful search direction, and dp mz * n and dpmax are 
the minimum and maximum dp z *, respectively. 

Searches then exhaust bin(i) before considering 
6in(i -f 1). 

3.1 Completeness 

Unfortunately, this path planning method is not 
complete, i.e., it cannot guarantee finding a solution 
even if one exists. Though certainly undesirable, this 
lack of completeness does not seem unreasonable since 
researchers have thus far been unable to develop algo- 
rithms which achieve both completeness and practical- 
ity for reasonably difficult yet practical path planning 
problems for more than a few degrees of freedom. We 
sacrificed completeness in exchange for the possibil- 
ity of solving some practical yet potentially difficult 
problems as quickly as possible. 

3.2 String Tightening 

Once a safe path is found, it may be modified to 
reduce the joint space trajectory length of the path. 
This process is referred to as string tightening [7]. 
Since the path planner produces discretized paths, the 
objective during string tightening is to reduce the fol- 
lowing cost function: 




iV-l 


t=i 




E (*,.(.■+ u-tyi))* 
j = 1 


(3) 


where: 

= the joint space trajectory length 
N = number of knot points in path 
n = number of dof 
0j(i) = knot point for joint j 
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The tightening algorithm involves examining each 

”“br,S' 0 ;7r ,on! p ' oduc “ ,he mos? 

1 Make no changes to the knot points. 

2. Modify the second knot point so that the three 
knot points are straight in the robot’s joint space 
(if not already so). P 

chel^n I e fhi bilily fi° f ° ption ; 2 must be determined bv 
caladZtm configuration for interference. These lo- 

DathuntiZ "ft contlnued al °ng the length of the 
from ™ b ' ob “ i "' d 

4 Implementation 

p-le I rn^ ddltl0n t0 havin s been implemented for sin- 
the C I e * r a Vk ^ , ro ^ 0 ^ P at ^ planning problems for 
scribed in fhi ed ' \ be path P lannin S strategy de- 
ASALoLth ni paper ha f , been implemented for the 
AbAL path planning problem described in Section 2 

The programs are written in C and utilize sect on" of 

by , Schi .T I*!- They .Uo i„Xm«h. 

The nolvf d devei °P ed b - v Hamlin and Kelley [9, 101. 
caused t y n ®P!. representatlon scheme was chosen be- 

5 lT accurate modeling of the robots and 

atfvelv ?f aC e V n the ' vorkceI1 w Hile enabling rel- 
atively fast interference checking. Paths are visually 

usTs 242 sea S rr n h g f im ? ation [H] The implementation 

S rrilS io d n 1M ' trans ‘" d 5 bi “ «■' **'•=" 

tended Ce t« thiS i Wa f a . Preliminary implementation in- 

p1a^ner for V t a be a \ e 9i h T e pos * ibI t usefulness of the path 

pl * n " i " s ‘ >tobiem - s °'"' 

• Nodes were not modeled. 


• In the ASAL, panels are installed (the first set af- 
ter the 60 strut). These panels were not mod- 
study) ?t f ° r ° ne partlcular strut as a case 

5 Results 

The path planner quickly found paths for the first 

?h at ,m« ,n n ‘ h f" “ litlk '" ,ssib '' mterference’at 
that stage Due to symmetry, the assembly of the 

remaining 81 struts can be accomplished using only 21 
unique trajectories for the Merlin with the appropriate 
carnage and turntable positions for each strut P The 

? 02 h s P a " ner f ' Vas , able t0 > d feasible paths for all 
102 struts with solution times ranging from 1 to 30 
minutes on a SparcStation 1. with the majority 

of solution times in the 2 to 5 minute range. J y 

n ..? e , 6lS . strut ls P os siblv the most difficult from a 
path planning perspective due to the confined location 
of the goal position and due to the presence of an 
installed panel above the goal position. Although this 
implementation generally ignored the panels /panel 

oHh/nf 6 ^ “ a u ° bst c de f ° r this strut - In s Pite 
• , * P anek a P atb was found without requiring any 
intermediate carriage/turntable positions The path 
found for this strut is illustrated in Figure 6 P 

mentaU^n follow^ COmments re S ardin S this imple- 

• The path planner has no trouble with goal posi- 
tions placing the load or robot in very close prox- 
unity to obstacles. 

* Inmk ath ^jP lan " er Performs well even with a large 

"Z e 'rfk 0bstacl l s ,- For f xam Ple, the final few 
struts of the assembly involve over 100 workspace 
obstacles. The additional collision checks P re- 
quired near the end of the assembly seem to 
increase execution time by a factor of approxi- 
mately two. 


33 






• The paths found typically include segments which 
are obstacle boundary tracing. Because of the 
close tolerances involved, it is not practical to sim- 
ply model the objects larger than actual size to 
provide a safety margin since doing so would often 
result m an unsol vable problem. 


• Panels and nodes were not modeled. As a result, 
some of the paths might collide with the panels or 

ui the paths were used in an actuaJ assem- 

bly- T ” lS COuld be remedie . <1 simply by modeling 
the panels and nodes and including them in the 
collision checking routine. Due to the small size 
of the nodes it is expected that including them 
would have little impact on the difficulty of the 
path planning problems. Although the panels 
will typically represent significant obstacles to be 
avoided, a strut for which the panels would seem 
to interfere the most was solved with the relevant 
panel modeled. 

• In a few cases the path planner was not able to 
solve the problem quickly in the forward direction 
but could quickly solve the problem in the oppo- 
site direction. Although a very confined goal po- 
sition makes it likely that solving in reverse mav 
prove easier, trial and error was the only sure way 
to decide which direction would vield better per- 
formance. 


• Return paths for the robot after inserting a strut 
were not planned. 


6 Conclusions and Future Work 

This implementation of the path planner for the 
ASAL assembly task illustrates the potential useful- 
rroccrf p 1 at * 1 planning technique developed at 
i-oT 5 .7 P?. , vin ? the Practical and potentially verv 
aimcult ASAL path planning problem. Based on the 
results of this study, additional work appears war- 
ranted towards applying this planning technique to 
the path planning problems at the ASAL. Some par- 
ticular issues which would need to be addressed before 
paths created by the planner could be executed on the 
actual hardware are as follows: 

• Nodes and panels need to be modeled. 

• An improved string tightening algorithm or an al- 
ternate method of path modification is required 
to provide paths with adequate clearances. This 
could be done by modifying the current string 
tightening cost function to include a penaltv 
on clearance or by utilizing the path found bv 
the planner as input to a potential fields based 
smoothing algorithm. 
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ABSTRACT 

This paper describes the implementation of a dis- 
turbance rejection controller for a 6-DOF PUMA ma- 
nipulator mounted on a 3-DOF platform. A control 
algorithm is designed to track the desired position and 
attitude of the end-effector in inertial space, subject 
to unknown disturbances in the platform axes. Exper- 
imental results are presented for step, sinusoidal, and 
random disturbances in the platform rotational axis 
and in the neighborhood of kinematic singularities. 

Robotic manipulators have been proposed as a 
means of reducing the amount of extra vehicular activ- 
ity time required for space station assembly and main- 
tenance. The proposed scenario involves a robotic ma- 
nipulator attached to some mobile platform, such as a 
spacecraft, satellite, or the space station itself. 

Disturbances in the platform position and attitude 
may prevent the manipulator from successfully com- 
pleting the task. This work explores the possibility 
of using the manipulator to compensate for platform 
disturbances. 

The problem of controlling a robotic manipulator 
on a mobile platform has received considerable atten- 
tion in the past few years. Joshi and Desrochers [1] 
designed a nonlinear feedback control law to carry out 
tasks (with respect to the robot base frame) in the 
presence of roll, pitch and yaw disturbances in the plat- 
form axes. Dubowsky, Vance, and Torres [2] proposed 
a time-optimal planning algorithm for a robotic manip- 
ulator mounted on a spacecraft, subject to saturation 
limits in the attitude control reaction jets. Papadopou- 
los and Dubowsky [3] developed a general framework 
for analyzing the control of free-floating space manipu- 
lator systems. Most recently, Torres and Dubowsky [4] 
have presented a technique called the enhanced distur- 
!i? nC ir map * man ipulator trajectories that reduce 
the effect of disturbances in the spacecraft position and 
attitude. 

One common assumption in the literature is that the 
disturbance signal is exactly known. If this is the case, 
then the end-effector location can be calculated with- 
out relying on direct end-point sensing. However, this 
assumption is invalid if there is a significant delay in 
the platform position and attitude measurements, or if 
the kinematics of the platform are not well known, or if 
the platform is a non-rigid structure. In the more likely 
case that only the nominal platform location and up- 
per bound on the disturbance signal are known, direct 


end-point sensing is needed to measure the end-effector 
location. 

1 The Jacobian and Singularities 

The inverse Jacobian is ill-defined for certain manip- 
ulator configurations. This section presents an alter- 
native mapping, called the approximate pseudoinverse 
Jacobian, which is defined for all manipulator configu- 
rations. 

The Jacobian maps differential changes in joint po- 
sition to differential changes in Cartesian position and 
orientation according to the following relationship: 

du = J(q)dq (1) 

where du E 9? 6 is the differential Cartesian displace- 
ment vector (linear and angular), q E is the vector 
of joint positions, dq E 3? n is the vector of differential 
joint displacements, and J E S 6xn is the Jacobian ma- 
trix. 

For the PUMA, the Jacobian matrix is simplest 
when expressed in frame 6: 


V 3l9 = 


(d$ + de)C^e 

(d$ + de)Ss6 

** 5^5 + 06^56 + ^ 75*56 

£56 

C56 

0 


di + cl$Sg d? 

&6 + < 15^6 a 6 

0 0 

0 0 

0 0 

1 1 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

—S 7 

C 7 Ss 

- 1 

0 

—Cb 

0 

c 7 

S 7 Ss 


( 2 ) 


The following compact notation will be used to denote 
the matrix 6 J 3| $: 


6 


^3,9 = 


B 0 ' 
D E 


(3) 


where £, D y and E are 3 x 3 submatrices of the Jaco- 
bian. 

The inverse Jacobian, when it exists, can also be 
written in block matrix form 
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5 J — 1 — 
J 3,9 ~ 


B 

D 


0 

E 


M- 


B" 1 0 1 

E- l DB- x E- 1 J 

. ( 4 ) 

The PUMA has three singularities. The first is re- 
ferred to as the Arm Fully Stretched singularity. This 
singularity occurs whenever the arm switches between 
the Ilex and the noflex configurations. 

The second singularity corresponds to the Hand 
Over Head singularity. The Hand Over Head configu- 
ration corresponds to changing between the right and 
left configurations. 

The third singularity is the IVrisf singularity , and 
occurs when the arm switches between the flip and 
noflip configurations. 


2 Approximate Pseudoinverse Jaco- 
bian 

The usual method of dealing with singularities of 
the Jacobian is to avoid them. This approach is not 
applicable to the disturbance rejection problem since a 
sufficiently large disturb ance could force the manipu- 
lator into a singular configuration. Also, the manipu- 
lator must avoid not just singular points , but singular 
regions , since the norm of J -1 becomes very large in 
the neighborhood of a singularity. 

The pseudoinverse Jacobian is often used to over- 
come the difficulty of J being a nonsquare matrix, and 
is defined as 


j 


{ 






m < n 
m = n 
m > n 


( 5 ) 


Clearly, this method of computing J * does not ad- 
dress the issue of singularities since it still relies on 
matrix inversion. A more general approach to com- 
puting the pseudoinverse uses singular value decompo- 
sition. This has one serious drawback, which is the 
high cost of computing the singular value decomposi- 
tion (SVD). The SVD algorithm uses a series of House- 
holder transformations to reduce the input matrix to 
diagonal form. Since this is an 0(N 3 ) operation, find- 
ing the SVD for the 6x6 Jacobian matrix can be too 
costly to implement in real-time. The alternative pre- 
sented in this section is called the approximate pseu- 
doinverse Jacobian, and is denoted by J *. 

The basic idea behind the approximate pseudoin- 
verse is to use the partitioned form of J and perform 
the SVD on the submatrices B and E. This reduces the 
number of computations by a factor of four, since two 
3x3 singular value decompositions is an 0(2(N/2 ) 3 ) 
operation. 

The definition of the approximate pseudoinverse Ja- 
cobian is 

Jt = [ - - jjtDBt ] ( 6 ) 

where B , D, and E are defined as in (3). 

Several properties of the approximate pseudoinverse 
are stated below. 


1. J* = J -1 when J is nonsingular. 

2. J* does not satisfy the Moore-Penrose conditions 
when J is singular. 

3. Let dp, d<f> G ft 3 be the linear and angular com- 
ponents of du, respectively, and let dq\,dq 2 E 3J 3 
be the components of dq. Then, the approximate 
pseudoinverse solution is 

[&H-IU *][*] < 7 > 

If J is singular, the approximate pseudoinverse 
finds the minimum norm solution as if dp and 
d<f> were decoupled ; that is, dq = J*du minimizes 
\\Bdqi - dp|| 2 and || Edq 2 - d^|| 2 . 

3 Behavior Near Singularities 

Figure 1 compares the 2-norm, or the maximum sin- 
gular value, of Jt (solid curve), J * (dashed curve), and 
J -1 (dotted curve) in the vicinity of the Hand Over 
Head singularity. 


1 

8 



Time (tec) 



Figure 1: 2-Norms of (solid curve), J * (dashed 
curve), and J -1 (dotted curve) Near Hand Over Head 
Singularity 

The discontinuities in \\J*\\ 2 and ||J*|| 2 occur when 
the smallest nonzero singular value, cr r , falls below the 
threshold value, Setting <r m i n to a relatively 

small value will shrink the width of the “well” about 
the singular point, thus extending the range over which 
jt = j- 1 and J* = J -1 . The side-effect is that the 
norm will be very large and highly discontinuous near 
the singularity. By the same token, setting (Twin to 
a relatively large value will reduce the discontinuity 
in the norm by increasing the width of the singular 
region. A threshold value of <r mtn = 0.1 was used to 
generate Figure 1. 
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4 Bound on Approximation Error 

The pseudoinverse and approximate pseudoinverse 
Jacobians are identical only when J is nonsingular. In 
order to characterize the difference in behavior at a 
singularity, some measure of the approximation error 
is needed. A reasonable way to measure the approxi- 
mation error is to see “how close” J* is to being a true 
generalized inverse using the following norm: 




Consider the matrix 


(8) 

Table 1: Computation Times for an d 

Subtracting J yields *^1 ,e 


JJ'J = 


BB'B 0 

D-(I-EE')D(I-B'B) EE'E 


Coordinate Frame 
k 

Computation Time | 

*7-1 

j 3.b 







i 

1.31 ms 

25.31 ms 

6.38 ms 

2 

1.19 ms 

25.3l ms 

6.25 ms 

3 

MEEim 



4 

KMYM7XW 

24.65 ms 

5.98 ms 

5 


24.65 ms 

5.98 ms 

6 




7 

iimza 



8 


Kft&vnn 


9 


24.65 ms 

5.85 ms 

E 





JJ'J- J= 


I 

0 

7 


0 

I -EE' 

B'B 0 
0 7 


B 0 
D E 


( 9 ) 


When both B and E are singular, the approximation 
error is bounded as follows: 




i 

0 


0 

7- EE ' 


r B 

0 ‘ 

‘ 

D 

E 



I -B'B 
0 


0 

7 


( 10 ) 


< imi 2 

If B is nonsingular, a less conservative upper bound 
can be found: 


||jj»j-7|| 2 


7 0 

0 I -EE' 


< ll^lla 


0 0 

0 E 

J 2 

(ii) 


Likewise, when E is nonsingular the upper bound re- 
duces to 


-•'i, = 

< 


fB O' 

1 - B'B 0 1 

[ 0 0 

0 1 J 


H L J L ~ * JII2 

l|S|| 2 (12) 


Finally, if both B and E are nonsingular, the approx- 
imate pseudoinverse is identical to the pseudoinverse: 


l J</ * 7 - / llj = ° (13) 


5 Computation Time 

Table 1 compares the computation times of the the 
inverse, pseudoinverse, and approximate pseudoinverse 
Jacobians for each coordinate frame. As predicted, the 
approximate pseudoinverse is about four times faster 
to compute than the pseudoinverse. 

The inverse, pseudoinverse, and approximate pseu- 
doinverse Jacobian solutions were implemented in the 
C programming language using the GNU 1 gcc Version 
2.2.2 compiler. The data in Table 1 was collected by 
timing the software on a Motorola MVME 147SA-2 
Single Board Computer. 

6 A Kinematic Control Law for Distur- 
bance Rejection 

Consider a 6-DOF PUMA manipulator mounted on 
a 3-DOF platform. The goal is to maintain the desired 
position and attitude of the end-effector with respect 
to the inertial reference frame (frame 0), subject to 
arbitrary disturbances in the platform axes. The fol- 
lowing information is assumed to be available: 

1. 9 E 8i 6 , the PUMA joint positions 

2. rj 0 E Sfi 3 , the nominal platform joint positions 

3. 6 E 3ft 3 , the maximum deviations from the nomi- 
nal platform joint positions 

4. °uo,« € Si 6 , the inertial end-effector location 

Two factors contribute to the motion of the end- 
effector: the differential displacement of the PUMA 
joints, which can be measured, and the differential dis- 
placement of the platform joints, which is unknown. 
Let 6 denote the disturbance signal and let dv be the 
component of the end-effector motion caused by the 

Copyright (C) 1989, 1991 Free Software Foundation, 
Inc., 675 Mass Ave, Cambridge, MA. 



differential displacement of the platform joints. Then, 
the differential end-effector displacement can be writ- 
ten as 

°dtio >E = °*/3, £?(*?<> "H Q)dO + dv 

= 3 R(rj 0 + 6) 3 Ja, E (9)d$ + dv (14) 

Note that coordinate frame transformations have been 
applied to isolate the dependence of the PUMA Jaco- 
bian on the platform joint positions. 

A discrete-time model of the system will now be 
derived by approximating the differential quantities 
in (14) with displacements. The underlying assump- 
tion here is that tie sampling period, AT, is suffi- 
ciently mall (Le., the sampling rate is much higher 
than the bandwidth of the system). 

Define Au* as A where the sub- 
script k denotes the kth sample step. In the limit as 
AT goes to icto, the displacement Ati* equals the dif- 
ferential du: 

lim Aujfc = du (15) 

Similarly, A 0* — ► dO and A Vk — * dv as AT — ► 0. 
Therefore, the discrete-time approximation is 

du « Au* = ut — tijb_i 

dO « A0* — 9 k — i 

dv « At;jk sr Vk — Vk~i (16) 

and the discrete version of (14) is 

°«t- °«i_i = ImVo + StfJzAWMk + Avk (17) 

where the subscripts denoting the reference and veloc- 
ity frames of du have been dropped to avoid confusion 
with the time index. 

Let u tid be the desired position and orientation of 
the end-effector along some specified trajectory. The 
control objective is to drive the end-effector to this 
position and orientation: 

°u* — ► °ti<* as k — ► oo (18) 

Ideally, the control objective could be achieved in min- 
imum time by computing the PUMA joint displace- 
ments A Oa needed to cancel out the inertial-space er- 
ror. However, exact cancellation would require com- 
plete knowledge of the disturbance signal. The next 
best solution then is to compute a A 0& which approx- 
imately cancels out the inertial-space error. With this 
goal in mind, the proposed control law is 

= 3 Jl E (0k)ln.(v>)K e (°U4 - °u k ) (19) 

where K c E 3ft 6 x 6 is a matrix of control gains. Equa- 
tion (19) will be referred to as the J* control law in 
the sequel. 


A simple expression for the closed-loop system can 
be derived by assuming that there is a one period delay 
in the control actuation: 

A0* +1 = AO* (20) 

°«* - %_! = §R(»j 0 + 6 k ) 3 J 3iB (0 k ) 

( 21 ) 

3 4, B (^-i)?R(»?o)^( 0 Ud - V-i) + At* 

In order to simplify this expression, define the quantity 

M M _j £ 3R(»?o + ^) 3 J3 lE (^) 3 4 B ( (, ‘-l)oR-(»7o)^ 

( 22 ) 

Rewriting (21) in terms of Mk t k-u it is easy to see 
that the closed-loop system is linear with time- varying 
coefficients; 

V = (J - M k , t - i)°u t _i + M m _ x 0 u d + At* (23) 

A block diagram of the closed-loop system is shown in 
Figure 2. 



Figure 2: Block Diagram of Closed- Loop System 

6,1 Design Parameters 

The selection of the control gain is greatly simpli- 
fied by restricting K c to be a scalar times the identity 
matrix: 

K e = Jb c 7, 0 < k e < 2 (24) 

The parameter k c controls the spectral radius of Aft. 
For example, if k c — 0.5, then the eigenvalues of M* 
will lie on a circle of radius 0.5 in the A-plane (or at 
zero, if J is singular). 

It is straightforward to choose a stable k c if 6 is 
known a priori . (Recall that b is the vector of maxi- 
mum deviations in the platform joint positions.) Let A 
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denote the spectrum of the matrix 3 R (r ) 0 + $) 

By invoking the slowly time- varying condition, a can 
be approximated as follows: 

a « sup arg(A) (25) 


and k c is calculated as 


c Vtan 2 a + 1 


(26) 


The selection of a m i n is essentially a trade-off be- 
tween tracking accuracy and the norm of the control 
signal. The selection of cr min should be based on the 
desired upper bound on the norm of A 6 d) which in 
turn is dictated by the saturation limits of the joint- 
level controller. 




7 Experimental Results 

Three sets of experiments focused on the time re- 
sponse of the closed-loop system for step disturbances 
in the platform joints, sinusoidal disturbances in the 
platform joints, and random disturbances in the plat- 
form joints. Here we present only the results for the 
step disturbance. 

This section analyzes the time response of the 
closed-loop system for 10° and 30° step disturbances 
in the platform rotation. For each case, the control 
gain K c was set to identity. 

7.1 10° Step Disturbance 

Figure 3 shows the inertial-space errors errors when 
a 10° step disturbance is applied to the platform rota- 
tional joint. The linear (X, Y , and Z) components of 
the error are shown in the upper plot and the orienta- 
tion error in the lower plot. The components of A 64 , 
the control vector, are plotted in Figure 4. 


Figure 3: Position Error ( X - solid curve; Y - dashed 
curve; Z - dotted curve) and Orientation Error for 10° 
Step Disturbance in Platform Rotation 



Figure 4: Control Signals (A0d(l)i A^(4) - solid 

curves; A0 d (2), A 0 d (5) - dashed curves; A0<*(3), 
A0d(6) - dotted curves) for 10° Step Disturbance in 


Platform Rotation 



Maximum Overshoot 

4% Settling Time 

X 

1.527 xl0 +u cm 

1.54 s 

Y 

3.825 xl0 +u cm 

0.84 s 

Z 

6.366 xlO -1 cm 

1.70 s 

<t>* 

3.503 xl0 +u deg 

1.00 s 


Table 2: Maximum Overshoot and 4% Settling Time 
for 10° Step Disturbance in Platform Rotation 

Table 2 lists the maximum overshoot and 4% set- 
tling time for the X, Y, Z, and orientation errors. The 
4% settling time refers to the time required for the er- 
ror to enter and remain within ±€ of zero, where e is 
4% of the peak absolute error. 

7.2 30° Step Disturbance 

The inertial-space errors and control signals for the 
30° case are shown in Figures 5 and 6. The maxi- 
mum overshoot and settling time for each coordinate 
are displayed in Table 3. 



Figure 5: Position Error ( X - solid curve; Y - dashed 
curve; Z - dotted curve) and Orientation Error for 30° 
Step Disturbance in Platform Rotation 






Control Signal (deg) 10*0) 




6: Control Signals (A0 d (l), A0d(4) - solid 
A0 d (2), A0 d (5) - dashed curves; A0 d (3), 
A0<j(6) - dotted curves) for 30° Step Disturbance in 


Figure 

curves; 


Platform Rotation 


8 Behavior Near Singularities 

Figure 7 shows the vector of open-loop control sig- 
nals near the Arm Fully Stretched singularity. The 
minimum singular value parameter, <r min , was set to 

0.1. At this value of the control in the direc- 
tion of the workspace boundary becomes very weak 
approximately 30° from the singular point. This pre- 
vents the end-effector from getting too close to the 
workspace boundary. Consequently, the manipulator 
will not switch between the flex and noflex configu- 
rations while the J * controller is running. 

If the parameter <r min is sufficiently small, however, 
the width of the singular region will be reduced to 
the point where the control signal for joint 6 (A0<*(3)) 
could drive the arm through the singularity. This 
may lead to an undesirable “chattering” behavior, in 
which the arm rapidly oscillates between the flex and 
noflex configurations. 

9 Summary 

Several important conclusions can be drawn from 
the experimental results. 

1. The relative stability of the closed-loop system 
is a function of the amplitude of the disturbance 
signal. 

2. The relative performance of the controller is a 
function of the frequency of the disturbance sig- 
nal. 



Maximum Overshoot j 

4% Settling Time 

X 

1.737 xl0 +1 cm 

1.97 s 

Y 

1.706 xl0 +1 cm 

2.43 s 

Z 

1.253 xl0 +1 cm 

1.66 s 

4>e 

2.055 xl0 +1 deg 

2.08 s 


Table 3: Maximum Overshoot and 4% Settling Time 
for 30° Step Disturbance in Platform Rotation 



Figure 7: Behavior of l/det(J) and Open-Loop Con- 
trol Signals (A0<*(1), A0<*(4) - solid curves; A0<*(2), 
A0<*(5) - dashed curves; A0<*(3), A0<*(6) - dotted 
curves) Near Arm Fully Stretched Singularity 


In other words, the J * controller is like a high-pass 
filter; the lowest frequency components of the distur- 
bance signal are attenuated the most. 

3. The control in certain directions becomes very 
weak near singularities. 

This implies that there may be an unavoidable tracking 
error in the “forbidden” directions when the arm is at 
or near a singularity. 
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Abstract 

This paper presents the results of controlling A PUMA 560 Robotic Manipulator 
and the NASA shuttle Remote Manipulator System (RMS) using a Command Gener- 
ator Tracker (CGT) based Model Reference Adaptive Controller (DMRAC). Initially, 
the DMRAC algorithm was run in simulation using a detailed dynamic model of the 
PUMA 560. The algorithm was tuned on the simulation and then used to control 
the manipulator using minimum jerk trajectories as the desired reference inputs. The 
ability to track a trajectory in the presence of load changes was also investigated in 
the simulation. Satisfactory performance was achieved in both simulation and on the 
actual robot. The obtained responses showed that the algorithm was robust in the 
presence of sudden load changes. Because these results indicate that the DMRAC 
algorithm can indeed be successfully applied to the control of robotic manipulators, 
additional testing was performed to validate the applicability of DMRAC to simulated 
dynamics of the shuttle RMS. 


1 Introduction 

Direct adaptive control offers the potential for uniform control of robotic manipulators 
in the presence of uncertain flexibilities, changing dynamics due to unknown and varying 
payloads, and nonlinear joint interactions without explicit parameter identification. 

One such direct adaptive algorithm that is especially attractive for robotic control is 
the direct model reference adaptive controller (DMRAC) discussed in [1-3]. This adaptive 
algorithm is very appealing for robotic control because of the following features: 

• asymptotically zero output error with all states bounded, 

• lack of dependence on plant parameter estimates, 

• direct applicability to multiple input*multiple output plants, 

* ThU P a P €r “ ba » cd u P<>n research performed under NASA grant NAGW-1333 and under NSF grant 
ECS-9 111565. 
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• sufficiency conditions which are independent of plant dimension, 

• control calculation which does not require adaptive observers or the need for full 
state feedback, 


• ease of implementation, and 

• successful experimental validation. 


This procedure has been previously used to control a single link flexible robotic joint 
aiid a nonlinear model of a two link Puma [4,5]. In view of the excellent tracking results 
demonstrated in these papers, it was concluded that this adaptive algorithm should be 
used to control an actual Puma arm. This effort has consisted of two mai n thrusts: namely, 
control of the representative simulation model developed in [6], and the transition of the 
tuned algorithm to the actual robotic arm. 

Because results indicated that the performance of the DMRAC algorithm was robust 
with respect to representative load variations, additional applicational studies were initi- 
ated using the NASA shuttle Remote System (RMS). 


2 Direct MRAC Development 

2.1 Basic algorithm 

The linear time invariant model reference adaptive control problem is considered for a 
plant which is described by 

MO = A p x p {t) + B p u,(t) ( 1 ) 

MO = C p x p (t) (2) 

where x p (t) is the (n x 1) plant state vector, Up(t) is the (m x 1) control vector, y p (t) is 
the ( q x 1) plant output vector, and A p , B p are matrices with appropriate dimensions. 

The objective is to find, without explicit knowledge of A p , and B p , the control u p (t) 
such that the plant output vector y p (t) approximates “reasonably well” the output of the 
following (and usually lower order) reference model: 

*m( 0 = A m x m (t) + B m u m (t) (3) 

ym(0 = C m X m {t) (4) 

The MRAC algorithm is given by [lj: 

MO = JMOliMO - MO] + x.(0*»(0 + K u {t)u m (t) ( 5 ) 

with the gains K e (t) t K z (t), and K u (t) being adaptive. The adaptive gains are concate- 
nated into the matrix K r (t ) which is given by 

K r (t) = [fr e (t),tf,(0.tfu(0] (6) 
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and the vector r(t) is defined by 


r (<) = [j/m(0 -y P (0,arm(0. U m(0] T (7) 

Then 

u p (t) = K r (t)r(t) (8) 

The adaptive gains are obtained as a combination of an integral gain and a proportional 
gain as shown below: 

Kr(t) = K p (t) + Kf{t) (9) 

*p( 0 = [y«(0 - y P {t)\r T (t)f (io) 

Kj{t) = [ym(t)-y p {t)]r T {t)T ( 11 ) 

Sufficient conditions for stability derived for a constant model command in [2]. These 
conditions require that the matrices T and T be respectively chosen as positive definite 
and positive semidefinite, and that the plant be almost strictly positive real (ASPR), 
that is, for the plant represented by the triple (A p , B P ,C P ) there exists a matrix K t (not 
needed for implementation) such that the fictitious stabilized plant described by the triple 
( A p — B p K e C p , B p ,C p ) is positive real. If these sufficient conditions hold, then all states 
and gains are bounded and the output error vanishes asymptotically. 

The adaptive control of plants that are not ASPR is a more difficult problem when 
utilizing the CGT based MRAC laws. BarKana [3] suggested augmenting the plant with 
parallel dynamics such that the augmented plant is ASPR in which case the previously 
described adaptive controller may be utilized. To illustrate this concept, consider the 
non-ASPR plant described by the transfer matrix 

G P {s) = C p (sl - A p y l B p (12) 

Then, choose a matrix H(s) such that the augmented plant transfer matrix 

G a (s)=H- l (s) + G p (s) ( 13 ) 

is ASPR. In [3] it is shown that G a {s) will be ASPR provided that 

• H(s) itself is ASPR 

• H(s) stabilizes the closed loop output feedback system with transfer function. 

[/ + G p (s)^(s)]- 1 G p (s). 

An easily implementable version of H (s) which has had extensive use is 

H(s) = if(l + s/so) 

resulting in: 

G,(«) = -. + ?. m +G p (s) (14) 
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where D = K~ l . 

Unfortunately, the error which is ensured to be stable is not the true difference between 
the plant and the model. Rather, the error is now the difference between the outputs of 
the augmented plant and the model. However, in [3] it is shown that if the plant is output 
stabilizable via high gain output feedback, then || D || may be chosen to be small. Thus, 
the augmented plant error will be approximately equal to the original plant error. 

One procedure for eliminating this output error is to also incorporate the supplemen- 
tary feedforward of (14) into the reference model output as shown in [2]. To illustrate this, 
denote the augmented plant output as z p where 


z p — Vp + Hu p 

(15) 

= y p + H[K x x m 4 - K u u m K e e v \ 

(16) 

and 


H = D/{ 1 + sr). 

(17) 

In a similar manner, define an augmented model output as 


z m Vm 4" H [K z X m + R" u U m ] 

(18) 


Now, for adaptive control of the augmented plant, consider the error between the aug- 
merited plant and model outputs. Thus, 


{^m Zp) — !/m yp HK € e M 

J/m yp HK e (z m Zp) 

or z m -z p = {I + HK € )~ l {y m -y p ) (19) 

Consequently if as in [2] z v is forced to follow z m , then y p will follow y m . 

3 Puma Model Development 

In order to test the performance of the Direct Model Reference Adaptive Controller (DM- 
RAC), an accurate non-linear coupled model of the PUMA manipulator was needed. A 
full explicit dynamic model of the PUMA 560 manipulator, derived by Armstrong, Khatib, 
and Burdick [6], was selected. The formulation of the PUMA model was computationally 
efficient using 25% fewer calculations than a 6 degree of freedom recursive Newton-Euler 
method. The algebraic formulation of the model also allowed for the easy addition of a 
load by modifying the link 6 mass, center of mass, and inertia parameters. 

Figure 1 shows the six rotational joint axis, { z \ , . . . , Zq}, for the PUMA 560. Only the 
rotational, z lf axis are shown in the figure. Positive rotations follow the right hand rule - 
counter-clockwise looking down the z axis. The six joint of the PUMA 560 are as follows: 

• Joint 1. A vertical rotation about the base, z\. 
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• Joint 2. A horizontal rotation about the shoulder, z 2 . 

• Joint S. A horizontal rotation about the elbow, z 3 . 

• Joint 4 • A twist of the wrist, z 4 . 

• Joint 5. A inclination of the wrist, z 5 . 

• Joint 6. A twist of the mounting flange, z$. 

The position of the manipulator in Figure 1 illustrates the zero position. Note that when 
Joint 5 is at zero, axis z 4 and z 6 coincide. 

The dynamic equations of motion used to model the PUMA are: 

A{$)6 + B{9)\99\ + C{0)[9 2 } + g{9) = T (20) 

where 

A(6) is the 6x6 positive definite kinetic energy matrix, 

B(9) is the 6 x 15 matrix of coriolis torques, 

C(6) is the 6x6 matrix of centrifugal torques, 

g{6) is the 6 vector of gravity torques, 

9 is the 6 vector of joint accelerations, 

[ 99 \ is the 15 vector of velocity products, where 

[Ms, Ms,..., Ms, Ms,..., 

Me.Me ] 7 

[9 2 \ is the 6 vector of squared velocities, where 
and T is the 6 vector of joint torques. 

The above model can be cast into state space form by solving Equation (20) for 9 and 
choosing the following 12 x 1 state vector, 

x T = [ 9 T y \ (21) 

where 

9 = 
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Thus, the state space model is as follows, 

0 = v 

The controlled output vector for the plant was 


Vplant = 0 + av (22) 

where, a is a diagonal 6x6 matrix of velocity weighting factors. 

This velocity term is present to help remove high frequency oscillations caused by the 
controller. The maximum allowable torques (in n-m) were [97.6, 186.4, 89.4, 24.2, 20.1, 


4 Implementation Issues 


4.1 Reference Model 

The first decision to be made in implementing the DMRAC algorithm is the choice of 
reference model order. If one chooses the order too low, then excessively large gains may 
occur even in a well-tuned controller. This may produce greater than desired accelerations 
in the robot arm joints resulting in joint torque saturations leading to poor model following. 
If one chooses the order too high then excessive response delays may be incurred. For the 
PUMA 560, an independent second order reference model was selected for each of the six 
joints. This is not unreasonable since in a PUMA 560, as with many manipulators, the 
mass matrix is approximately diagonal for all 0 making the system nearly decoupled. 
Thus, for each joint, the reference model transfer function was: 


where 


y«,(a)/u m ,.(s) = 


tv: 


s 2 + 2&tu n< .s + w 2 . 


Wn. = 5 

and 

6 = 1 * = { 1 , 2 , 3 , 4 , 5 , 6 }. 

Critical damping was selected so as to reduce the possibility of joint angle overshoots. 
This conforms to a standard safety feature of robot arm controllers which tends to avoid 
obstacle collisions. Of course, once the choice of critical damping is made, the choice of 
natural frequency governs the speed of model response to inputs. A choice of w n = 5 
yields a 90% rise time of about 0.8 sec. 
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4.2 Command Generation 


For Testing purposes a minimum jerk trajectory was generated through the following 
positions at the noted times. 



Joint 

Position (deg) 



1 

2 

3 

4 

5 

6 

Time (sec) 

0 

-45 

180 

0 

45 

90 

0 

90 

-90 

90 

45 

0 

45 

6 

0 

0 

180 

0 

90 

90 

13 

0 

-45 

180 

0 

45 

90 

18 


The resulting angular paths for each joint were then used as the reference model com- 
mands u m> .(i). 

4.3 Bias Introduction 

For the PUMA 560 manipulator, the origin of the coordinate system should be such that 
the adaptation gains have a non-zero excitation throughout the range of interest. For 
example, assume that in order to maintain an output of y p = [0, ...,0] T , a non-zero input, 
Up, is required. However a zero command to the reference model, u m = [0, ...,0] T , will 
result in a zero model output and a zero state vector. Thus in this case e v = y m — y p 
will also be zero, and the vector, r(t), defined by (3.10) will be zero resulting in a zero 
control. Since the plant requires a non-zero control to maintain a zero output, the DMRAC 
algorithm requires a non-zero error signal in order to apply a non-zero control which will 
result in a steady-state error at the zero output position. 

If the reference model coordinates are shifted by a constant bias term, then a zero 
command to the reference model, u m = [0, ...,0] T will produce non-zero outputs for the 
model state and output vectors which, in turn, will produce a non-zero command to the 
plant. This bias term is subtracted from the model command, u m , and the plant output, 
y p , as follows, 

u m (t) = u m (t) - qua, (23) 

!/ P (t) = &(*) ~ ( 24 ) 

where u m (t) is the original model command in the original coordinate system, u m (t) is 
the new biased model command to be applied to the model dynamics, y p {t) is the actual 
plant output, y p (t) is the new biased plant output to be used to form the error signal, 
and qua, is a constant bias term. For robotic manipulators, qua* Has units of radians 
and should be selected such that a new plant output of y p = [0, ..., 0] correspond to an 
equilibrium position. By examining the zero position of the robot, Figure 1, it is clear 
that y p = {0, 0, 0, 0, 0, 0} is not an equilibrium. However bias of, 

9km = {0, 90, 90, 0, 0, 0} degrees (25) 

will shift the zero position to that shown in Figure 2. 
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( 26 ) 


4.4 Feedforward Design 

The feed-forward filter dynamics for joint t is given as, 


A (5) 


1 + TiS 


where K d . is the DC gain and r, is the time constant. 


5 PUMA Simulation Results 

Ja this section, we briefly discuss the tuning process and present plots of a simultaneous, 
six joint response of the PUMA 560 under DMRAC control. 


5.1 Tuning 

Once the reference model has been chosen, one must choose values for the various DMRAC 
parameters. Specifically these are 

T = proportional gain weighting matrix, eq. (10) 

T = integral gain weighting matrix, eq. (11) 

D = plant/model feedforward gain, eq. (17) 
r = plant/model feedforward time constant, eq. (17) 
a = 6 vector of plant rate feedforward gains, eq. (22) 

For the fully centralized DMRAC algorithm with the plant derivative output term and 
the supplementary feed-forward in the reference model and plant, there are 1182 param- 
eters to be selected as shown in Table 5. At first, this number seems very intimidating, 
but as will be shown, the number of tuning parameters can be greatly reduced by some 
simplifications and by adjusting the parameters in groups rather than individually. 

The most drastic reduction in the number of tuning parameters can be achieved by 
forcing the integral and proportional adaptation weighting matrices, T and T to be diag- 
onal. This reduces the number of tuning parameters from 1182 to 78. 

The reference model dynamics have 12 tuning parameters, six w m 'a and six f.’s. 

It is customary in robotic applications to tune controllers for critical damping so that 
there is no overshoot. Overshoot may cause a robot end effector to penetrate the surface 
of the work environment. The undamped natural frequency terms, w n . are chosen such 
that the reference model will have a specified step response. Typically, the reference 
model dynamics are chosen such that they are “reasonable” for the plant to follow since 
the DMRAC algorithm will try to force the plant to follow the model output. For the case 
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of a PUMA 560 Manipulator, all of the w ni were initially set to 5.0. The model’s dynamic 
parameters can be changed as needed if the robot is having problems tracking the model. 

Initially, the plant output derivative weights, a, are set to zero. These weights are used 
to remove high frequency components from the plant control signal, u m , and should only 
be used when needed as they will' affect the transient response. 

The feed-forward filter has 12 tuning parameters, six gains K 4 . and six time constants 
r,. A good first choice for the is approximately one-tenth the model time constant. 

Initially the value of r, were all set to 0.1 s, and the six DC filter gains were set to 
1.0. Increasing the filter gain was seen to typically improve the tracking performance. 

The diagonal components of T and T were initially, all set to unity. A reasonable 
method of tuning a DMRAC controller is to start the plant at an equilibrium position and 
apply small step inputs. After a reasonable performance is achieved with the step inputs, 
the DMRAC should be fine tuned using typical plant trajectories. 

If the closed loop system is very sensitive to initial conditions, start with small steps 
as described above, let the system reach steady-state, and then save all of the DMRAC 
controller state information (integral adaptation matrix, K\\ reference model state vector, 
x m ; and the filter state vector) to be used as initial conditions for the next run. This will 
significantly cut down the adaptation time required for the gains to reach their steady-state 
values. 

In order to compare the tuning results, some criterion must be established. For this 
excessive, the goal was to kept the peak model following errors small and to keep the error 
trajectory as close to zero as possible. Small errors were tolerable during motion. It was 
also desired to achieve zero error in steady-state. 

The step response with the initial tuning values was sluggish for Joints 1, 4, 5, and 6 
with overshoot and oscillations. Joints 2 and 3 settled into their steay-state values quickly 
but with very large steady-state errors. The process used to complete the tuning was as 
follows: 

1. Refine the tuning for a 10 degree step from the equilibrium position. 

2. Using the refined parameter values, move the robot to the shutdown position of 
Figure 12 and save the DMRAC internal state values at that position for use as 
initial conditions. 

3. Refine the tuning for a 10 degree step from the shutdown position using the initial 
conditions from Step 2. 

4. Refine the tuning from typical min-jerk trajectories from the shutdown position. 

The final tuning parameter values after Step 4 are shown in Table 6. The weighting 
matrix values for Joints 1, 2, and 3 differ from the weighting matrix values for the last three 
joints by a factor of about 100 which reflects the mass/inertia difference between the upper 
arm and the wrist. The weighting matrix values which are multiplied by the “Xm,” products 
are about a factor of seven lower than the values multiplying the u Xmi n products since the 



second state variable of each decoupled reference model had a higher peak value in a 
transient. The Joint 1, 2, and 3 reference models have am undamped natural frequency 
of 4.0 rad/sec where the wrist model used 7.0 rad/sec which again reflected the inertia 
difference between the upper arm and the wrist. The feed-forward filter values were set 
to Kd = 6.0 and r = 0.1 for all joints. The alpha values were increased from the initial 
values of zero to damp out some high frequency oscillations. 


5.2 Response 

Initially the PUMA arm is in the [0 0 0 0 0 0] position. The final shutdown position was 
(0, -45, 180, 0, 45, 90) degrees as shown in Fig. Stimulation results of the PUMA 560 
dynamics responding to the tuned DMRAC controller are displayed in Figure 4. 

Note that the model following is excellent for all 6 joints. Furthermore it was observed 
that all joint torques were smooth and below their saturation limits. In addition, for 
this specific case, the use of the feedforward component did not significiantly affect the 
response, although in other cases (eg. step response) use of the feedforward resulted in 
significant improvements. 


6 PUMA Experimental Results 

Because the simulator results of the previous section indicated that DMRAC should be 
useful for robot control, a set of experiments was performed on an actual PUMA 560 
manipulator. The tuning process was similar to that described in the previous section. 
All parameters were initalized at those values from the simulation studies. Only minor 
variations were required. Final values are in Table 1. 

Examples presented illustrate performance of the DMRAC for tracking various trajec- 
tories in the presence of static and dynamic load changes. In all cases the robot starts at 
the shutdown position and follows a trajectory which finishes at the shutdown position. 

6.1 Three Joint Trajectory Tracking Study 

The trajectory listed in Table 1 is very similar to the one used in the simulation (Section 
5). The arm first moves to a straight up position, curls up, and then moves back to the 
safe position. The wrist joints remain locked in their shutdown positions of {0.0,45.0.90.0} 
degrees. 

The response to the first trajectory is shown in Figure 5. The response is quite good. 
The effects of stiction can be seen on Joint 2 at t = 15 seconds in Figure 5. Figures 6-8 
show the model following error and the link torques for Joints 1, 2, and 3 respectively. 
Figure 6b shows that the Joint 1 torque signal was quite noisy. This noise did not have 
a physically detectable effect on the actual arm motion. Typically one can feel or hear a 
noisy torque signal on the actual arm. 
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The stiction effect mentioned above for Joint 2 can also be seen in figure 4a at t = 15 sec 
near the ‘X’ at the peak error location. When stiction grabs a joint, the error ramps up 
as does the torque (Figure 4b). F v 


6.2 Static Load Changes 


Thu. section descnbes the ability of the DMRAC algorithm to adjust to static load vari- 
a ■ons The trajectory of Table 3 will be run with different loads in the gripper. The 

algorithm will first be allowed to adjust to the load, and then the trajectory will be 
started. 


The wrist joints remained locked in their shutdown positions of {0.0, 45.0, 90.0} de- 
grees. Five different loads were run for the trajectory - 0kg, lkg, 2kg, 3kg, and 4kg. 

Figures 9 and 10 show the response for Joints 2 and 3 respectively. The numbers on 
the plots are to help identify which curve represents which payload. For Joint 3, the peak 
errors vary from 2.4390 degrees for the no load case to 3.9972 degrees for the 4 kg load 
case. The load changes make up only about 50% of the error. The other 50% is due to 
the adaptation to the changing arm dynamics. For Joint 2, the peak errors are around 
0.8-1.0 degrees. As with Joint 3, the portion of the error due to the load change for Joint 
2 is small compared to the no load case. 

,._ For J ° int j the error signals did not vary by more than 0.1 degrees between the five 
different load cases. 


6.3 Dynamic Load Changes 

To illustrate the effects of dynamic load change, the trajectory of Table 4 was considered. 
While running the same trajectory, various loads were added to the gripper while the robot 
were in motion. The same loads used in the previous section were employed. The wrist 
joints remained locked in their shutdown positions of {0.0, 45.0, 90.0} degrees. Note: The 

lkg wid 4 kg loads were added at about t = 6.76 seconds and the 2 kg and 3kg loads were 
added at about t = 7.34sec. 

Figure n shows the model following error for Joint 2 for all loads. The numbers on the 
graphs indicate which peaks in the error plots match up with the various loads. This figure 

shows that the DMRAC algorithm has a good load disturbance rejection. The transient 
period only lasts about 2 seconds. 

F^ure 12 shows the error for Joint 3 for the various loads. Joint 3 suffers more 
W :‘\ a A ° ad dls J urbanc e having a peak error of almost 5 degrees when the 4 kg load is 
a ed. Again, the transient period is roughly 2 seconds. After the transient, good tracking 
performance was achieved with the additional loads. 

As with the static load case, the model following errors for Joint 1 did not vary by 
more than 0.1 degrees. 
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7 Adaptive Control of the Shuttle RMS 

7.1 Introduction 

Because of the previous demonstrated capabilities of DMRAC, consideration was given 
to its application to the NASA shuttle Remote Manipulator System (RMS). This system 
experiences damped oscillations of the end effector after the motion control input by the 
shuttle operator has been removed [7], It is desired to design a controller that will take 
control of the RMS, after the operator releases the motion control joystick, and increase 
the damping of the oscillations. Linear models have been developed for three manipulator 
orientations expected to be encountered during normal payload handling [7]. This section 
discusses work on a direct model reference adaptive controller design for the RMS based 
on the adaptive algorithm discussed in Section 2. 

7.2 Linear Plant and Feedforward 

The three linear plants are 3-input, 3-output with 6 states. The plants all have a feedfor- 
ward compensator H(s ), since they are not ASPR. Three types of algorithm feedforward 
were examined in the course of this work, 
static: 

tf(s) = diag3{dll,d22,d33} 

1st order: 

H(s) = diag3{dll/ (r„s + 1)} 

2nd order: 

H(s) = diag3{d rt /[(r lri s + 1) * (r 2l ,s -I- 1)]} 

The scalar feedforward provided the best results (based on work to date) and was used 
for all presented simulations. It was found that, for the scalar feedforward, the combined 
plants (plant 1,2 or 3 in parallel with feedforward) were all ASPR for: 

0.125 < d < 1.0 

where H(s) was diag3 {d,d,d}. That is, the closed loop system formed from the inverse 
of H(s) in negative feedback with the respective plant, had all the characteristic roots in 
the left half plane for: 

1.0 < d~ l < 8.0 

After simulations were performed with many of the possible combinations of values within 
this range, it was found that d = 0.25 provided the best results for all three orientations. 

No ASPR analysis was performed for the 1st and 2nd .order feedforwards and the po- 
sition 1 plant. The stability of the adaptive algorithm for the dynamic feedforwards was 
somewhat a function of the adaptive gains, T and T, for given feedforward time con- 
stants and gains. Tuning for the dynamic feedforwards was difficult, and very little to 
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no improvement over the uncontrolled system could be achieved. The dynamic feedfor- 
ward compensator might possibly provide improved results with further work on the time 
constants and gains. 


7.3 Reference Model 

Originally, a reference model were developed for each position using a LQR design based 
on the uncontrolled plant at the respective orientations. The model used the uncontrolled 
plant’s B, C, and D matrices with the model A matrix formed as follows: 

•^ m i (-^p, -6p, * K%) 

where if, is from the LQR design for the i-th orientation. 

Satisfactory control of the plants could not be achieved by the adaptive algorithm using 
these models. A new reference model was then developed using two dominant eigenvalues 
from the original LQR model for position 1. The new model has a damping ratio of c and 
a natural frequency of 1.0 r/sec so that: 

hrn( ,s ) = + 2 $S + l) 

This new model was utilized as the reference for each plant output, that is: 

H m {s) = d\a.g3{h m (s),h m (s),h m (s)} 

^m(^) * U m (s) 


7 .4 Simulation Sequence 


The sequence for simulation represented a 3.0 second perturbation followed by use of the 
controller to dampen out oscillations. 

The three plant outputs were: 


Yl = shoulder yaw 
Y 2 = shoulder pitch 
Y 3 = elbow pitch 

The plant inputs were limited to 0.7 deg/sec.. 

In order to simulate the perturbation, the following control 
uncontrolled plant: 


sequence was input to the 


ui = 0.7 0 < t < 1.5 
t*i = -0.7 1.5 < t < 3.0 
u 2 = u s = 0. 

The resulting plant states at the end of this perturbation were the plant initial conditions 

° r * . er ® imulatlons - Th e plant outputs at the end of the perturbation were the 

model initial condition (s). The model rate initial condition (s) were set to zero. 
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^•5 Simulation Results 

Figures 13 a,b,c show the position 1 outputs for the following parameters: 

T = diag9(6000, 10, 6000, 1 , 1 , 1 , 1 , 1 , 1) 

f = diag9(l, .000001, 1, 100, 100, 100, 100, 100, 100) 

D = diag3(.25, .25, .25) 
model damping = 0.15 


Note that all controlled outputs decay faster than do the outputs with no control 

Figures 14 a,b,c show the position 2 outputs using the above position 2 controller 
unmg parameter and Figures 15 a,b,c show the position 1 outputs using controller 
parameters tuned for position 2. For these cases the differences between the controlled 
ana uncontrolled responses were not remarkable. 

“ d ° thCr eXperiments show that a satisfactory level of control can be 
achieved by the MRAC with tuning tailored for the individual positions. Attempts to 
develop one set of controller tuning parameters that would provide satisfactory control for 
all three positions were not successful. 


8 Conclusions and Recommendations 

h summary, the DMRAC algorithm has been found to be an effective robotic control 
algorithm in both simulation and on the actual robotic manipulator. Its performance was 
robust with respect to static and dynamic load variations and also disturbances. 

At present the DMRAC is being considered for all six joints of the actual PUMA and 
turther tuning with dynamic feedforward is begin considered for the shuttle RMS. 
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Table 1: Parameter Values for 3 Joint Trajectory Tracking 


f 

a n 
c z 

20 

40 

40 

(diag 

» 

140 

20 

200 

component) 


30 

200 

30 


U u w 
u m 

140 

200 

200 

T 

c * 

30 

30 

40 

(diag 

x m 

200 

30 

400 

component) 


60 

400 

60 



200 

400 

400 

joint 

1 

2 

3 

Model 


10 

10 

10 


<r 

1 

1 

1 

Feed 


6 

6 

6 

Forward 

T 

0.05 

0.05 

0.05 

alpha 

a 

0.02 

0.02 

0.02 


Table 2: First Three Joint Tracking Test Trajectory 


Knot 

Point 

Joint Positions (deg) 

Time 

1 

2 

3 

(sec) 

0 

0 

-45 

180 

- 

1 

-90 

-90 

90 

6 

2 

0 

0 

180 

8 

3 

0 

-45 

180 

6 


Table 3: Static Load Change Trajectory 


Knot 

Point 

Joint Positions (deg) 

Time 

1 

2 

3 

(sec) 

0 


-45 

180 

- 

1 

0 

-45 

180 

3 

2 

45 

0 

0 

10 

3 

0 

-45 

180 

10 


Table 4: Dynamic Load Change Trajectory 


Knot 

Point 

Joint Positions (deg) 

Time 

1 

2 

3 

(sec) 

0 

— 

-45 

180 

- 

1 

0 

-45 

180 

3 

2 

45 

-90 

90 

10 

3 

0 

-45 

180 

10 


6l 


Runs 







Parameters | Description II Values 


24 x 24 integral weighting matrix 576 


24 x 24 proportional weighting matrix 576 


Undamped natural frequency for Joint i model 6 


Damping ratio for Joint i model 6 


6x6 diagonal plant derivative weighting matrix 6 

DC gain of Joint t supplementary feed-forward block 6 

Time constant of Joint i supplementary feed-forward block 6 


Total 1182 





u " 140 160 110 


T 

(diag 

component) 




Joint 




6 6 6 6 6 6 


0.1 0.1 0.1 0.1 


a I 0.0.35 0.02 0.02 0.01 0.01 0.01 


Figure 1: PUMA 560 Coordinate Frame 
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Figure 2: PUMA 560 in Stable Equilibrium 



Figure S: Shutdown Position 






turn, mc 


Figure 4: Simulation results for all 6 joints 



Figure 5: Plant and Model Output for First Trajectory, 
(a) Joint 1. (b) Joint 2. (c) Joint 3. 
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(a) Joint 1. Model Following Error 



Figure 6: Joint 1 Data for First Trajectory, (a) Model following error, 
(b) Joint Troque. 



Figure 7: Joint 2 Data for First Trajectory, (a) Model Following error, (b) 
Joint Torque. 
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Figure 8: Joint 3 Data for First Trajectory, (a) Model following error, (b) 

Joint Torque. 



Figure 9: Joint 3 Static Load Model Following Error 
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Figure 10: Joint 3 Static Load Model Following Error 


JotM 2, Model Following Em m 



Figure 11: Joint 2 Dynamic Load Model Following Errors 
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Joint 3. Model Following Errors 
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Figure 12 


Joint 3 Dynamic Load Model Following Errors 





TIME (sec) 



TIME (sec) 

Figure 13: Shuttle RMS - Orientation 1 outputs for MRAC tuned to this 

position 
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Figure 15: Shuttle RMS - Orientation 1 outputs using the MRAC tuned for 
orientation 2 
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1 Introduction 


The Explorer Platform is a Modular Mission Spacecraft; it has several sub- 
units that are designed to be replaced on orbit. The Goddard Space Flight 
Center Robotics Lab undertook an experiment to evaluate various robotic 
approaches to replacing one of the units; a large (approximately* 1 meter by 
1 meter by 0.5 meter) power box. The hardware (see figure 1) consists of 
a Robotics Research Corporation K-1607 (RRC) manipulator mounted on 
a large gantry robot, a Kraft handcont roller for teleoperation of the RRC\ 
a Lightweight Servicing Tool (LST) mounted on the RRC, and an Explorer 
Platform mo ckup (EP) with a removable box (MMS) that has fixtures that 
mate with the LST. Sensors include a wrist wrench sensor on the RRC*, and 
Capaciflectors [Vranish91] mounted on the LST and the MMS. There are 
also several cameras, but no machine vision is used. The control system lor 
the RRC is entirely written by Goddard [Leake91]; it consists of Ada code 
on three Mulitbus I 386/387 CPU boards doing the real-time robot control, 
and C on a 386 PC processing Capaciflector data. The gantry is not moved 
during this experiment. 

The task is the exchange of the MMS; it is removed and replaced. This 
involves four basic steps: mating the LST to the MMS, demating the MMS 
from the EP, mating the MMS to the EP, demating the LST from the MMS 
Each of the mating steps must be preceeded by an alignment to bring the 
mechanical fixtures within their capture range. 

Two basic approaches to alignment are explored; teleoperation with the 
operator viewing thru cameras, and Capaciflector based autonomy. To eval- 


2 ALIGNMENT 


uate the two alignment approaches, we ran several runs with each approach, 
and recorded the final pose. Comparing this to the ideal alignment pose 
gives accuracy and repeatability data. In addition, the wrenches exerted 
during the mating tasks were recorded; this gives information on how the 
alignment step affects the mating step. 

There are also two approaches to mating; teleoperation, and impedance- 
based autonomy. The wrench data taken during mating using these two 
approaches is used to evaluate them. 

Section 2 describes the alignment results, Section 3 describes the mating 
results, and finally Section 4 gives some conclusions. 


2 Alignment 


The two alignment tasks are aligning the LST for mating with the MMS, 
and aligning the MMS for mating with the EP. Two methods were used for 
each task; teleoperation, and Capaciflector-based autonomy. 

For teleoperation, we used the Langley rate control algorithm. One experi- 
enced operator performed all the runs. The Kraft hand- controller acts like 
a 6 DOF joystick; the rate of the RRC tool frame is proportional to the 
displacement of the Kraft from a reference frame. In a traditional joystick, 
there is a centering spring force returning the joystick to the reference frame; 
in the Langley algorithm, this centering force has a constant magnitude, not 
proportional to the displacement. This allows wrench feedback to be added 
to the centering force without operator confusion. On the RRC, there is 
a Cartesian impedance algorithm using the wrench sensor, that, makes the 
RRC tool frame behave like a pure damper; it relaxes when any force is 
applied. Thus if the tool is against a surface, and the operator pushes the 
hand- controller into the surface, the hand- controller commands a constant 
rate, which is turned into a constant force by the damper algorithm on the 
RRC. At the same time, the wrench sensed by the RRC wrench sensor is fed 
back to the motors on the hand-controller. Wrench feedback ratios ol 1:1 
can be achieved with this algorithm. Since the alignment task is primarily 
free-space positioning, the wrench feedback was low for this task, to mask 
noise and errors in the gravity model of the loads carried by the wrench 
sensor. Three cameras were used during teleop; one on the LST (only used 
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2 ALIGNMENT 


for mating the LST to the MMS), one giving an overall view of the RRC 
and EP , and one giving a good view of the MMS mounted on the EP. 

For Capaciflector-based autonomy, there are two Capaciflectors mounted on 
the LST, and six on the MMS. Alignment is a 6 DOF task; this is easily 
accomplished with the six sensors on the MMS. For the LST, there is no 
way to place six sensors to get a full alignment. So a sequential approach 
is used; first the two sensors are leveled against a surface, then they find an 
edge, then they find a bump along the edge, etc. Each step in the sequence 
can find two degrees of freedom; we used a 7 step sequence to help eliminate 
errors. 


For each task, a perfect” goal pose was defined manually (using rulers and 
direct vision to align the LST and MMS). Then 10 runs lor each combination 
of task and approach were made, recording the final pose for each run. The 
accuracy is defined as the mean error between the “perfect” goal pose and 
the 10 actual poses; the repeatability is the standard deviation ol the same 
error. For the MMS, there are not enough visual cues to allow the operator 
to align the MMS with the EP without contact. So the MMS. TELEOP task 
did both alignment and initial contact; the MMS -CAPACIFLECTOR task 
did not contact. The following table summarizes the results, and the scatter 
plots in figures 2 thru 5 show the raw data. 


run 

time (sec) 

translatio 

accuracy 

n (mm) 
repeat 

rotation ( 
accuracy 

radians) 

repeat 

LST.TELEOP 

53 +- 10 

7.30 

5.29 

0.02789 

0.07483 

LST.CAPACIFLECTOR 

123 +- 0.5 

5.69 

0.45 

0.01883 

0.00096 

MMS.TELEOP 

42 +- 10 

7.49 

15.5 

0.01920 

0.04911 

MMS.CAPACIFLECTOR 

61 +- 20 

2.88 

11.7 

0.01092 

0.03554 


For both tasks, teleoperation is significantly faster. The LST task shows 
a 23% improvement in translation accuracy from teleop to Capaciflector. 
and a 32% improvement in rotation accuracy. The repeatability is signif- 
icantly better with the Capaciflector; a factor of 11 for translation. 77 for 
rotation. For the MMS task, the accuracy shows a 61% improvement in 
translation, and 43% improvement in rotation. The repeatability shows a 
24% improvement in translation, and 27% in rotation. Remember that the 
teleoperation MMS task used contact for the final alignment, while the Ca- 
paciflector made no contact; the Capaciflector algorithm is more repeatable 
than the mechanical contact! 
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3 MATING 


For the MMS, it often took 1 or 2 incorrect contacts before the final correct 
contact was made. 

For the LST, the operator felt teleoperation was more reliable, while for the 
MMS, he felt the Capaciflector was more reliable. 


3 Mating 

Wrench data was logged for both LST and MMS mating, recording data 
from just before contact until after full contact. The fixtures in both tasks 
guarantee 6 DOF alignment. The LST is essentially rigid; the contact bp- 
tween the LST and the MMS fixture is basically a narrow cone, with a plate 
and two posts at the top for roll alignment. The clearance between the LST 
and the mating fixtures on the MMS is about 2 mm and 0.1 radians. The 
contact between the MMS and the EP is at three points ( before the screw 
is fastened); there is no clearance, but the contacts are actually spherical, 
so some misafignment is possible. No screws were tightened or latches fas- 
tened, to simplify data analysis. Two methods were used lor each mating 
task; C art esian impedance control, and teleoperation. The LST mating was 
then repeated using a more complex impedance control. For teleoperation, 
one experienced operator performed all the runs, starting irom the same 
starting point as the teleop aligment task. For autonomy, the start pose 
was representative of the final alignment pose using Capacifiectors. Thus 
the differences include differences in starting alignment as well as mating 
algorithm. 

For teleoperation, we again use the Langley rate algorithm. The wrench 
feedback and RRC damper gains are adjusted to give the best operator 
feel, while maintaining stability. The best operator leel is achieved when 
the joint stiffness is very low. If the joint stiffness is high, the Cartesian 
impedance loop has too work very hard to overcome it, and this shows up 
as instability at high wrench feedback gains. However, to use a low joint 
stiffness, we .must use gravity compensation torque, to keep the arm from 
sagging. Unfortunately, a design flaw in the RRC- analog servo hardware 
prevents us from using a torque command when carrying the MMS payload, 
so we had to use relatively high joint stiffness for the MMS. We repeated the 
LST runs using a softer joint servo with gravity compensation. The actual 


77 


3 MATING 


gains for each non are given in the appendix. 

The performance measures are the time from first contact to stable contact, 
rms wrench error, and the maximum wrench. The Z axis is the mating 
direction; wrench error in this direction is measured only after stable contact . 
The following table summarizes the results; figures 6 thru 13 show the raw- 
wrench logs for a representative run. 


run 

time 

TX 

ri 

TY 

ms wren 
TZ 

ich err' 
RX 

yr 

RY 

R.Z 

max w 
trail 

Tench 

rot 

LST-MATE.TELEOP 

14.6 

1.07 

3.84 

14.60 

0.95 

0.40 

0.12 

53.6 

2.48 

LST-MATE-AUTO 

6.0 

1.58 

1.58 

22.97 

1.65 

3.70 

0.14 

79.4 

8.29 

MMS -MATE-TELEOP 

32.2 

11.3 

30.7 

5 5. a 

4.69 

8.Q7 

2.38 

250.3 

27.8 

MMS .MATE-AUTO 

2.2 

r 3.82 

2. 38 

r n ri 

0.82 

1.77 

0.30 

82.1 

4.65 

LST-MATE-SOFT-TELEOP 

7.2 

4.49 

1.88 

6.14 

0.67 

1.26 

0.52 

51.6 

4.19 

LST-MATE-SOFT-AUTO 

4.5 

4 .fi() 

5.01 

15.7 

0.82 

1.77 

0.27 

105.4 

8.11 


For both tasks, teleop was slower than autonomy. 


For the LST task, teleop gave lower wrench errors, and lower m aximum 
wrenches. This is attributed to the fact that the operator used vision to 
refine the alignment as mating proceeded. 

For the MMS task, there is a factor of at least 5 improvement in wrench 
errors for autonomy; the off-insertion- axis portion can be attributed to the 
more accurate and reliable Capaciflector alignment, while the on- insertion- 
axis portion is due to the more stable mating algorithm. 

Note that the LST -MATE run was slightly unstable when fully mated, and 

that MMS -MATE-TELEOP typically made contact twice incorrectly before 

finally seating. 


The operator would always prefer to use the autonomous impedance algo- 
rithm, not teleop. The control system is actually more complex for teleop. 

The more complex gravity compensation control system gave faster times 
and lower on-insertion-axis forces, but higher ofi-msertion-axis wrench ei 
rors. It was more stable than the non-gravitv compensation system 
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4 CONCLUSIONS 


4 Conclusions 


For both tasks, the Capaciflector gave an improvement of at least. 25% in 
alignment accuracy and repeatability. For the LST task, the repeatability 
improved by a factor of at least 11. 

For the LST task, teleoperation alignment followed by teleoperation mat- 
ing gave lower wrench errors, by a factor of about 2. For the MMS task, 
Capaciflector alig nm ent followed by autonomous mating gave lower wrench 
errors by a factor of at least 5. These results are not conclusive; more work 
needs to be done to distinguish between the effects of initial alignment and 
mating algorithm. 

We anticipate significant reductions in the wrench errors and maximum 
wrenchs with future control system improvements, for both autonomy and 
teleoperation. 
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7 APPENDIX 


7 appendix 

LST -MATE -TELEOP uses teleoperation in the Langley mode, with the 
following parameters: 

Motion.Scale => (others -> 0.1), 

Wrench.Feedback *> (Active => TRUE, Scale => (others -> 0.25)), 
Joint.Servo *> 

(Joint.Servo.Label *> ANALOG.DAMPING.NOGRAV. 

Pos.Error.Action => CLIP, 

ADN.Vel.Gain => (6.0, 6.0, 6.0, 6.0, 6.0, 4.0, 4.0)), 

Cart. Impedance => (Active => TRUE, 

Bias «> (others => 0.0), 

Spring => (others => 0.0), 

Damper => (1000.0, 1000.0, 1000.0, 200.0, 200.0. 100.0)) 

LST.MATE-AUTO uses Cartesian impedance control, with the following 
parameters: 

Joint.Servo *> 

(Joint.Servo.Label => ANALOG.DAMPING.NOGRAV. 

Pos.Error.Action => CLIP, 

ADN.Vel.Gain => (6.0, 6.0, 6.0, 6.0, 6.0, 4.0, 4.0)), 

Cart .Impedance -> (Active => TRUE, 

Bias => (TZ => 40.0, others => 0.0), 

Spring ®> (others => 0.0), 

Damper => (4000.0, 4000.0, 4000.0, 1600.0. 1600.0. 400.0)) 

MMS-MATE.TELEOP uses teleoperation in the Langley mode, with the 
following parameters: 

Motion.Scale ■> (others => 0.05), 

Wrench.Feedback *> (Active => TRUE, Scale => (others => 0.125)). 
Joint.Servo => 

(Joint.Servo.Label *> ANALOG.DAMPING.NOGRAV, 

Pos.Error.Action => CLIP, 

ADN.Vel.Gain => (6.0, 6.0, 6.0, 6.0, 6.0, 4.0, 4.0)), 
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Cart .Impedance => (Active => TRUE, 

Bias ■> (others => 0.0), 

Spring => (others => 0.0), 

Damper => (4000.0, 4000.0, 4000.0, 800.0, 800.0, 800.0)), 

MMS _M ATE_A.UTO uses Cartesian impedance control, with the following 
parameters: 

Joint.Servo => 

(Joint.Servo.Label => ANALDG_DAMPING_NOGRAV . 

Pos_Error_Action => CLIP, 

ADN.Vel.Gain => (6.0, 6.0, 6.0, 6.0, 6.0, 4.0, 4.0)) r 
Cart .Impedance => (Active => TRUE, 

Bias => (TZ => 80.0, others => 0.0), 

Spring ■> (others => 0.0), 

Damper => (4000.0, 4000.0, 4000.0. 1600.0. 1600.0. 1600.0)) 

LST-MATE-SOFT-TELEOP uses teleoperation in the Langley mode, with 
the following parameters: 

Motion.Scale => (others => 0.2), 

Wrench.Feedback => (Active => TRUE, Scale => (others => 0.5)), 

Joint.Servo => 

( Joint.Servo.Label => PD.GRAV. 

Pos.Error.Action => CLIP, 

PDG.Stiffness => (9000.0, 3562.55, 2625.63. 2341.13. 341.75. 385.42. 80 
PDG.Damping = > (900.0, 580.59, 318.80, 281.04, 33.40, 37.04, 1 

Cart .Impedance => (Active => TRUE, 

Bias => (others =-> 0.0), 

Spring => (others => 0.0), 

Damper => (1000. 0, 1000.0, 1000.0, 100.0, 100.0, 100.0)). 


LST-MATE j>OFT-AUTO uses Cartesian impedance control, with the fol- 
lowing parameters: 

Joint.Servo => 

(Joint.Servo.Label => PD.GRAV, 


81 


7 APPENDIX 


Pos_Error_Action -=> CLIP, 

PDG.Stiffness => (9000.0, 3562.55, 2625.63, 2341.13, 341.75, 385.42, 80 
PDG.Damping => (900.0, 580.59, 318.80, 281.04. 33.40, 37.04. 1 

Cart .Impedance => (Active => TRUE, 

Bias => (TZ => 40.0, others =>0.0), 

Spring => (others =>0.0), 

Damper => (1000.0, 1000.0, 1000.0. 100.0. 100.0. 100.0)). 
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ABSTRACT 


A large-angle, flexible, multi-body, dynamic modeling capability was 
developed to help validate analytical simulations of the dynamic motion 
and control forces which occur while berthing of Space Station 
Freedom to the Shuttle Orbiter during early assembly flights. The 
paper describes the dynamics and control of the station, the attached 
Shuttle Remote Manipulator System, and the Orbiter during a berthing 
maneuver. Emphasis is placed on the modeling of the Shuttle Remote 
Manipulator Sytem in the multi-body simulation. The influence of the 
elastic behavior of the station and of the Remote Manipulator System 
on the attitude control of the station/Orbiter system during the 
maneuver is investigated. 
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STUDY PURPOSE AND PROCEDURE 


STUDY PURPOSE: Provide independent data to help JSC 
validate Space Station Freedom Program on-orbit assembly 
simulation analyses. 


STUDY PROCEDURE: Perform time simulations of the berthing 
of an Intermediate-build configuration of Space Station 
Freedom to the orbiter. 


A dynamic simulation of the berthing process is fairly complex since it 
involves the interaction of large, highly flexible components during a 
large motion maneuver while in orbit, where the components are 
subject to active control forces and gyroscopic, drag, and gravity 
gradient forces and moments. The complexities of the Space Station 
Freedom (SSF) assembly analytical simulator are such that it was 
advisable to develop independently a comparable tool to help validate 
the simulator. This paper is concerned with a description of a large- 
angle, multi-body, dynamic modeling capability developed to help 
validate the SSF program analytical berthing simulator which will be 
used to analyze each assembly flight. 

The berthing simulations are used to calculate the dynamic motion and 
control forces that occur while berthing early build configurations of SSF 
to the Orbiter during assembly flights when attitude control of the stack 
resides with the station control systems. The sixth assembly flight is 
the first flight that will use the station control systems rather than the 
Orbiter Digital Auto Pilot to maintain the attitude of the stack. Berthing 
during this sixth flight was selected as the validation simulation since 
the control systems of both the station and the Shuttle Remote 
Manipulator System (SRMS) are active during this maneuver. 
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MOTIVATION FOR SELECTION OF THE 
SIMULATION SCENARIOS 

Wanted to simulate a SSF assembly sequence 
which would exercise the: 

-SSF RCS attitude command system 

- SSF CMG attitude and momentum 
management system 

— Orbiter SRMS control system 

-Orbiter SRMS brake system 

-Orbiter SRMS flexible body 

- SSF flexible body 

-Environmental conditions during a typical orbit 

=> Selected with JSC the berthing of SSF- Stage 5 
to the Orbiter (MB-6 flight). 


The simulation scenarios were selected to capture the critical 
components of the space station assembly operation when attitude 
control of the stack resides with the SSF. The specific features which 
were to be exercised are outlined above. Since the MB-6 flight was the 
first flight to capture all of these features, it was chosen by JSC and 
LaRC to be the study scenario. 
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The figure shows the Space Station Freedom - Stage 5 configuration. 
The location of the avionics platform containing sensors which provide 
attitude and attitude rate information is indicated. The attitude can be 
controlled by firing jets, located on the top and bottom of the inboard 
station framework, at a constant force level of 25 lbs per jet, or by a set 
of four double-gimbaled CMGs, each with a capacity of 3500 ft-lb-sec, 
located on a platform close to the avionics platform. Also shown is the 
resource node, a pressurized shell attached to the station framework 
inboard of the alpha joint, to which a grapple fixture is mounted. 
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SSF/OR BITER APPROXIMATE DIMENSIONS 

FLIGHT MB-6: ORBITER BERTHED TO SC-5 STATION 
(UNITS -FEET) 



APPROXIMATE WEIGHTS 

SSF (SC-5) 145,000 lbs 

ORBITER 250,000 lbs 

(incl. payload) 

1000 lbs 


46.3 


122 


The f'gure shows the relative size and location of the stage 5 station 

a. S rb,ter ’ anc * tbe exten ded SRMS at the beginning of the simulation’ 
At this assembly stage, the SSF is over 150 feet in length. It has a 
weight of 145,000 lbs and the Orbiter, with the lab module in the carqo 
bay, has a weight of 250,000 lbs. The SRMS has a weight of only 
1 ,000 lbs. 1 
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TEA / BERTHING SIMULATION 


ESTABLISH TEA 


A 





s' :-v- 

X-,vVv'- 
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The scenario under investigation is the stage 5 assembly sequence 
depicted in the Figure. For the purposes of this study, this scenario is 
broken down into two simulations: 1) simulation of the Torque 
Equilibrium Attitude (TEA) maneuver and 2) simulation of the berthing 
of stage 5 to the orbiter. 

Before the first simulation begins, the photovoltaic (PV) arrays are 
feathered and the alpha joint is locked to minimize plume loads from the 
Orbiter jets during the final approach of the Orbiter before grappling 
occurs. The alpha joint remains locked during the entire berthing 
maneuver. The Orbiter approaches the station along the direction 
opposite the orbital velocity vector and flies in tandem with the station 
maintaining a distance of about 30 feet from the V guides in the cargo 
bay to the trunnion pins on the berthing adapter. The SRMS end 
effector grapples the station by snaring the grapple fixture located on 
the resource node. This is the start of the first simulation. The SRMS 
joint brakes are applied and the RCS jets are fired to move the station 
from a GG attitude to a computed Torque Equilibrium Attitude (TEA). 
Once TEA has been established, the second simulation begins. The 
brakes on the SRMS are released, the station RCS jets are inhibited 
Inn ®’ anc * att ' tuclG the stack is now maintained by the station 
♦k e D T?e5 ne ^ tum mana 9enient system. Berthing is accomplished by 
the SRMS joint motors which draw the station and the Orbiter together. 
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ESTABLISHMENT OF INITIAL TORQUE 
EQUILIBRIUM ATTITUDE (TEA) 




Stack at GG 


Stack at Initial TEA 


S hllHTnn?!l ibri ^ atti Vi de (TEA) is the average attitude which must 
ovjj lino nrh> • S0 that the net an 9 ular momentum accumulated 

orwta?otrnornn,> Sft p ^ esence . of Sravity gradient, aerodynamic and 
of tit SVJoscoP'c disturbances, is zero. The figure shows a schematic 
oHhe stack configuration at Gravity Gradient (GG) and at the initial 
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FIRST SIMULATION: 

ESTABLISH TORQUE EQUILIBRIUM ATTITUDE 

Use station RCS jets to move stack from attitude at 
grapple to a computed TEA. 

- Station RCS control system active, responding to 
commanded attitude. 

- Brakes with friction modelled are applied to all six 
SRMS joints. 

-Flexible representation of the SRMS long booms. 

-Flexible representation of SSF (36 mode model). 

-Orbital mechanics including aerodynamic and 
gravity gradient moments. 


The specific components exercised during the simulation of TEA 
maneuver are outlined here. The SSF Reaction Control System (RCS) 
is used to maneuver the stack to its TEA. During this maneuver, the 
SRMS brakes are applied. The SSF is modelled as a flexible body as 
are the upper and lower long booms of the SRMS. Aerodynamic and 
gravity gradients forces and torques are included. 




During the second simulation, the TEA changes since the inertia of the 

thfrPA^^f as , th ? Station is berthed t0 the Orbiter. The change in 
the TEA pitch angle dunng this simulation is shown in the Figure. 
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SECOND SIMULATION: 

USE SRMS TO BERTH STATION TO ORBITER 

- Station CMG control system active. 

-Flexible representation of the SRMS long booms. 
-Flexible representation of SSF (13 mode model). 

-Orbital mechanics including aerodynamic and 
gravity gradient moments. 

-Mass matrix and estimated TEA recomputed 
every second. 

-SRMS operational modes exercised: 

• Automatic mode 

• Manual Augmented mode 

• Position Hold mode 


The specific components exercised during the simulation of the 
berthing maneuver are outlined here. The SSF Control Moment Gyros 
(CMGs) are used to maintain the attitude of the Orbiter/SRMS/station 
stack. As in the TEA maneuver, the SSF and long booms of the 
SRMS are modelled as flexible bodies and the aerodynamic and 
gravity gradients forces and torques are included. The SRMS control 
system is active and is used to command the SRMS to berth the 
station to the Orbiter using Manual Augmented and Automatic Mode 
maneuvers. The changes in the TEA are estimated using a mass 
property estimator which computes the composite inertia of the stack 
as a function of SRMS end effector position and attitude. 
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SC-5 FE MODEL 

PV Arrays Feathered 


Number of Grid Points 

Number of Elements 

Number of Dynamic DOF 
after component mode reduction 

Number of Modes < 5 Hz 

Number of key modes selected 
for simulation 


= 36 (1st sim) 

13 (2nd aim) 


Jegr^r o f SSF °° nSiSted °< almost '5.000 

number of DO?warriduMd,o7m K°" 1 ? nt moda /educlion. the 
to 5Hz were ScSd M ° de shapes and '^enaes up 

™ e , ®P^ ce station structural dynamics were represented durina the TEA 
7 Hz To" c, y ose S ,o L^'Th’r^ Wfch ra "» ^ MS? 

were obtained for the model fixed at the grapple fixture point. 


107 



SSF ATTITUDE CONTROLLER DESIGN 



A simplified block diagram of the RCS and CMG control systems is 
shown. The attitude determination system (ADS), which measures the 
attitude and feeds this information back to the controller, is assumed to 
be accurate within the controller bandwidth so that a transfer function of 
unity is assumed for the ADS for the current simulations. The control 
system is designed for use in all configurations of the station covering a 
large range of inertias during the 3-year assembly process. To 
accommodate this wide range of system parameters, a mass estimator 
is provided to determine the on-orbit inertias and to adjust the control 
gains for acceptable performance. Normally the gains will not change 
significantly for a given flight configuration since the inertia matrix will 
remain nearly constant until the next assembly flight; however, during 
the berthing process, the system inertia matrix is continuously changing 
so that the gains are also continuously changing. A mass estimation 
program based on knowledge of the SRMS end effector location has 
been written to provide an updated system inertia matrix as berthing 
progresses. There are two bending filters (low pass filters designed to 
remove higher frequency components of the feedback position and rate 
signals) in the control design. 
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MOTOR SERVO CONTROLLER FOR EACH SRMS JOINT 



Each of the six SRMS joints is comprised of a reversible dc drive motor, 
a mechanical joint brake, an inductosyn tachometer, an epicyclic gear 
train and an electro-optical encoder and servo compensation as shown. 
The SRMS is telerobotically controlled from the aft cockpit of the Orbiter 
by way of translational and rotational hand controllers and control panel 
command inputs. Joint rate commands are sent from software 
algorithms resident in the Orbiter General Purpose Computer (GPC) to 
the joint servos by way of the Manipulator Control Interface Unit (MCIU) 
(not shown). The joint gear train applies the required torque to the 
SRMS/payload system. The encoders and tachometers provide 
measurement of the joint position and rate, respectively. 


no 




As shown, a digital joint rate command (in counts) is received from the 
GPC by way of the MCIU. This input rate demand is compared with the 
actual motor rate from a digital tachometer feedback to form an error 
signal. This error is then converted to an analog voltage signal and 
processed in through a trim integrator and a low pass filter. The 
purpose of the integrator is to provide a high gain at low frequencies 
needed to break motor and gear train stiction and to overcome small 
errors. The output of these analog electronics are summed with 
negative feedback of the analog tachometer signal. This continuous 
part of the tachometer is run through a high pass filter which serves to 
provide stability. In contrast, the purpose of the digital tachometer is to 
improve tracking accuracy. The analog voltage signal is then sent to 
the Motor Drive Amplifier (MDA) and a current limiter which results in 
attenuation of the voltage applied to the joint motor. The resulting 
motor rate is then passed to the gearbox to generate the required 
output torque. 

Nonlinear friction losses are modeled on both the motor or input side of 
the gearbox and on the output or joint side of the gearbox. The output 
friction models include both the joint friction and the gearbox friction. 


in 





JOINT DRIVE-TRAIN 



DRIVE 


Shown is a schematic of the SRMS joint drive-train. The optical 
encoder is physically mounted on the joint side of the gearbox and the 
tachometer is mounted on the motor side of the gearbox. The friction 
and freeplay which occurs in the drive-train have been modelled in the 
simulation. 
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NONLINEAR SERVO MODELS 


Gearbox Stiffness 

A ’'oeaitoo* 



t 


where, 

T a = gearbox transition torque (ft-lbs) 

Bj_ - gearbox transition angle (rad) 

Kq = gearbox linear stiffness (ft-lbs/ rad) 


Joint Friction/Stiction 



where, 

cSHerence in actual joint angle (deg) 
t c - coulomb torque level (ft-lbs) 
t s « stiction torque level (ft-lbs) 


"the old hairbrush model " 


The nonlinear gearbox model is represented by an asymptotic linear 
compliance and a quadratic stiffness relation at low torque levels. This 
stiffness is computed as a function of the backlash angle of the gearbox 
as shown. 


For the joint friction/stiction model, the friction torque is computed as a 
function of the joint angle as shown. When the actual joint angle is less 
than a stiction (static friction) reference angle, friction torque acts like a 
spring restraint. In this region, a steady torque produces a small 
rotational displacement. When an applied torque is removed, the joint 
returns to an equilibrium position. If the displacement exceeds the 
stiction reference angle, the friction torque drops to a coulomb torque 
level. While the joint rate remains positive, the friction torque is 
constant. If the rate becomes negative, the friction reference angle is 
reset and the model reverts to its spring-like behavior about the new 
reference point. This is the so-called “old hairbrush” model used in the 
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NONLINEAR SERVO MODELS (cont’d) 


Motor/brake Friction 



where, 

$ - motor rate (rad/sec) 
t c « coulomb torque level (ft-bs) 
t s - stiction torque level (ft-bs) 


The SRMS motor/brake friction/stiction is modelled as shown. When 
the joint is moving, the friction torque is equal to the coulomb (sliding) 
friction value. When the joint is not moving and the torque is larger than 
the stiction level, the friction torque is set to the stiction level. If the 
torque is less than the stiction level, just enough stiction is applied to 
the joint to make the net torque output zero. 


11U 



MUTLIBODY REPRESENTATION 


should «r pilch 

should* y wt 
swing-out (fixed) 





9 bodlos 

6 RMS joints modatod 
as sing# DOF revolul# joints 

111 -rigid 

- IbxlUa 


Nine bodies are used to model the complete multibody system as 
shown. The nine bodies include the orbiter, the seven links of the 
SRMS, and the SSF. Three components, the two long booms of the 
SRMS and SSF, are modeled as flexible bodies. Eight joints are 
defined to connect each of the bodies in the system. The swing-out 
joint at the base of the SRMS and the connection between the end- 
effector and payload are modeled as bracket (rigid) joints. The 
remaining six joints are modeled as single degree-of-freedom revolute 
joints, which accommodate the six degrees-of-freedom of the SRMS. 
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SIMULATION ARCHITECTURE 



The detailed architecture of the computational tool is shown. It consists of four 
major parts; multibody dynamics code DADS (Dynamic Analysis and Design 
System), the SRMS controller, the SSF ACS, and the MAIN program. The 
DADS code is used to generate equations of motion of the system, including the 
SRMS arm, the orbiter, and SSF. Each of these modules has its own 
integration routines and integrate its state equations at its respective integration 
step sizes. In order to synchronize the simulation process of different sampling 
rates, the MAIN routine was added to control the timing and program execution 
flow. 

For the SRMS controller, joint angles and rates from DADS, along with operator 
command inputs, are fed into the SRMS command algorithm to compute joint 
rate commands. The SRMS controller model calculates driving torques, based 
on the joint rate commands, which are then applied back to DADS through 
control elements. For the ACS, the DADS code provides the attitude and 
attitude rate of the stack to the ACS. Along with the commanded attitude, the 
ACS computes attitude errors and rate errors that are used to compute required 
commanded torques to be applied to the system. At the same time, the mass 
property estimator is used to estimate the inertia of the composite system. This 
estimated information is used to compute proper gain scheduling in the ACS 
and to update the commanded attitude. Depending on the type of actuator 
used, the commanded torques are converted to either RCS forces or to CMG 
torques which are fed back to DADS using control elements. Environmental 
disturbances from aerodynamics moments and gravity gradient torques, 
recomputed on the estimated inertia change, are also applied to the system 
through control elements. 
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SIMULATION RESULTS 


All sims performed using an SGI 4D/440 workstation 

ESTABLISH TORQUE EQUILIBRIUM ATTITUDE 
36 mode representation of SSF 
4 mode representation of SRMS 
TEA was established in 20 minutes real-time 

BERTHING 

13 mode representation of SSF 

4 mode respresentation of SRMS 

Berthing to within 2 feet was completed in 14 minutes 
real-time. 


For the TEA maneuver, the time integration was performed with a time 
step of 0.001 seconds and the computations took approximately 72 
hours of dedicated CPU time on an SGI 4D/440 workstation. The 
station structural dynamics were represented using 36 normal modes. 
A comparison was made with a simulation using ten normal modes to 
represent the station dynamics. The response differences were within 
the accuracy of the computation and thus the ten-mode model was 
deemed to be sufficient in representing the dynamics of the station. A 
conservative proportional damping level of 0.2 percent of critical 
damping was assumed for each mode. 

For the Berthing maneuver, the time integration was performed with a 
time step of 0.002 seconds and the station dynamics were represented 
using a 13 mode model. The maneuverwas completed in 14 minutes. 
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ATTITUDE CHANGE DURING TEA MANEUVER 



The TEA maneuver required an attitude change from a gravity-gradient 
position to an attitude orientation of pitch, yaw and roll of -22.1°, -7.4° 
and 3°, respectively. The resulting attitude-change time history is 
shown. The TEA of the stack was successfully established within 1 200 
seconds, a little more than a quarter of an hour. 


118 





The actual jet firing times for the upper and lower x-axis and z-axis jet 
clusters are shown. The RCS commanded torque is realized through 
the jet selection logic as firing pulses of approximately 0.2 to 1 second 
in duration. The RCS jets are inhibited from firing more often than once 
every 33 seconds to reduce structural dynamic response. 
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CHANGE IN POR MAGNITUDE DURING TEA MANEUVER 



0 500 1000 1500 2000 2500 

TMtfttq 


Shown is a plot of the relative magnitude of the position of the Point of 
Resolution (POR) (located at the SRMS end-effector) with the brakes 
engaged. The SRMS joints were slightly overloaded when the jets 
were fired to establish TEA and two of the wrist joints exhibited some 
brake slip during the first 750 seconds of the simulation. The largest 
position change computed during the jet firings was less than one inch 
and the resultant slip after the jet firings were completed was less than 
0.2 inches. The total position change is a combination of brake slip, 
arm flexibility and joint flexibility in the six-joint drive-trains. 
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POR LOCATION FOR BERTHING SCENARIO 

Orbiter Body Axis System (OBAS) Coordinates 



Results from the berthing simulation are shown. The SRMS POR is 
commanded to move along a three-point berthing trajectory. From the 
initial grapple point, position 1 the SRMS is commanded to move 4.5 
feet in the x and z-axis directions and 2 inches in the y-axis direction 
using an Operator Commanded Auto Sequence (OCAS) maneuver. 
From position 2, the SRMS is immediately commanded to move 25 feet 
vertically to position 3 using Manual Augmented Loaded mode z-axis 
translational hand controller inputs. Position 3, which is two feet from a 
full berthed position, is reached in approximately 800 sec. At this time, 
the Position Hold function is automatically enabled by the SRMS 
command algorithms to maintain its commanded position and attitude. 
Very little residual vibratory motion is observed following completion of 
the berthing maneuver. 


This simulation was conducted using translational and rotational end- 
effector rate limits of 0.14 ft/sec and 0.14 °/sec (coarse mode). 
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ATTITUDE CHANGE DURING BERTHING MANEUVER 



TMC (SEC) 


Cuter Aitfte Cuf u ntu te MteH, tm, Bo* 


The inertia of the stack, and thus the TEA, changes continuously during 
the berthing maneuver. The instantaneous mass of the stack is 
computed every second and used to modify the gains of the CMG 
momentum management control system. This information is also used 
to compute the current TEA which is subsequently applied to the CMG 
control system as an updated commanded attitude change. The 
change in attitude of the stack during the berthing maneuver is shown. 
The final Pitch, Yaw and Roll attitude of the stack is designated. 
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SIMULATION OBSERVATIONS 


ESTABLISH TORQU E EQUILIBRIUM ATTITUDE 
Minor SRMS brake slip occurs (continuous firing unacceptable). 
Influence of structural response on control Is small. 

Structural response and loads are small. 

RCS control functioned well. 

BERTHING 

CMGs saturate if maneuver starts with unbiased CMGs. 

SRMS Position Hold mode shows tendency toward Instability 
(SRMS upgrades not Implemented). 

Influence of structural response on control is negligible. 

Periodic update of stack Inertia Is required (trace of I matrix 
changes by as much as 10% during berthing). 


The following observations may be made about the simulated TEA and 
berthing maneuvers. The SSF RCS successfully established a TEA for 
the stack with only a small amount of SRMS brake slip. During the 
berthing maneuver, the CMGs were able to track the changing TEA 
while the orbiter and station were pulled together by the SRMS when 
given an initial bias. 

Although not shown, the Position Hold mode did exhibit a tendency for 
instability when left on for an extended period of time. This may be 
attributable to the known instability of Position Hold with massive 
payloads. An enhancement to Position Hold mode along with two other 
SRMS upgrades are presently being implemented by JSC to facilitate 
assembly operations. 

In both simulated maneuvers, the elastic behavior of the station and of 
the SRMS was found to have only a minor influence on the attitude 
control of the stack and the control loads caused only minor internal 
structural loads and structural response. 
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LaRC/MSFC Loosely Coupled Multibody 
Spacecraft Controls Research Facility 


Raymond C. Montgomery 
NASA Langley Research Center 
Hampton, Virginia 

and 

Dave Ghosh 

Lockheeed Engineering and Sciences Company 
Hampton, Virginia 
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LaRC/MSFC Loosly Coupled Multibody 
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PRESENTATION OUTLINE 

Problem 

Testbed Overview 

Preliminary Design, Analytical Models, and Simulation Tools 
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Space Station Berthing to the Space Shuttle 
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Multibody Controls Research Testbed 

A Joint MSFC/LaRC Testbed 
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Interim Research Testbed 
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MSFC Flat 
Floor Facility 





SSF Mobility Base Features 

Propulsion — 24 cold gas thrusters 

Sensors — 1 rate gyro, 2 accels, 2 laser scanners 

Attitude Control — 1 CMG 
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Preliminary Design (Continued) 
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Results from Preliminary Design 
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Cantilevered System Natural Frequenci 




















Photograph of the Testbed 

View from Mobility Base 
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Photograph of the Testbed 

View from Wall 
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Photograph Showing the CMG 
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Closed-Loop CMG Requirement 
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Y-Jet Firing Results for Conf. 

Y -Accelerometer 
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Wrist Joint Torque Input 

Gear Efficiency = 50% 
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Elbow Joint Torque Input 
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Status Report of RMS Active 
Damping Augmentation 


Mike Gilbert 

NASA Langley Research Center 
Hampton, Virginia 

and 

Martha E. Demeo 
VIGYAN 

Hampton, Virginia 
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STATUS REPORT OF RMS ACTIVE DAMPING 

AUGMENTATION 
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Workshop on Selected Topics in 
Robotics for Space Exploration 

March 17-18, 1993 
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Lockheed Engineering & Sciences Company / Hampton, VA 
Georgette W. Watts 
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ACTIVE DAMPING AUGMENTA 
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Joint LaRC / JSC program 




Reduced crew training time 
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Damping ratio improvement factor 



RMS ACTIVE DAMPING AUGMENTATION SCHEDULE 





SIMULATED PAYLOAD REDEFINITION 
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SENSOR AND ACTUATOR DEFINITION 





CSI CONTROLLER IMPLEMENTATION 
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MODIFICATION 




ADA CONTROL LAW DESIGN 
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Insure no adverse interaction with Shuttle FCS 
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STUDY STATUS 
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Phase 2b (completed September 22-24, 1 992) 

• Multi-point control law testing 

• Astronaut evaluation of RMS ADA 


ASTRONAUT EVALUATION 
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Unanimous recommendation to evaluate ADA with heavy payloads. 








THUOT / MT09241 2 / RUN 1 / MASC2 













ASTRONAUT RECOMMENDED SSF ASSEMBLY STUDY 
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INFORMATION NEEDED FOR FLIGHT DEMO DECISION 
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CONCLUDING REMARKS 
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- IBM for GPC software implementation of ADA 

- SPAR for accelerometer installation and overall approach 
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Autonomous Rendezvous and Docking 
for Space Station Freedom 
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James L. Garrison, Jr. and Stephen J. Katzberg 
NASA Langley Research Center 
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NASA Langley Research Center 


REQUIREMENTS FOR AR&D EXPERIENCE 
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automation and robotics, one element of which is a telerobotic 
research platform in space 

(Annual Report of the Langley GN&C Technical Committe - September, 1992) 
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Langley was asked to provide simulations and analysis of the mission, 
recommend algorithm enhacements, and perform ground-based hardware 
testing and simulation, if appropriate. 


COMET SPACECRAFT PERFORMANCE 
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Payload Volume: 15 cu. ft. 


AR&D MISSION PROFILE 



COMET 1 COMET 2 





ANALYSIS IN SUPPORT 
OF THE COMET AR&D MISSION 
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Definition of AR&D Standards 


ANALYTICAL MODELS IDENTIFIED FOR 
THREE REGIMES OF MISSION 
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MODEL FOR APPROACH/PROXIMITY OPS. 
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Sensor Noise 

Visibility of Laser Ranging and Machine Vision 
GPS Constellation Geometry 
Model of GPS Reciever 









MODEL FOR LOOSELY COUPLED CONFIGURATION 
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Enhancements: 

Sensor Model: Horizon sensor visibility, Magnetic Field Model 
Sensor Noise 

External (Random) Disturbances 
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MODEL FOR CONTACT DYNAMICS 
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APPLICATION TO SPACE INFRASTRUCTURE 
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1 1 . Participating in Commercial Space Ventures . Office of Commercial Programs, May 1 990. 

12. Tchoryk, P M Apley, D. J. f Jr., Conrad, D. J., and Dobbs, M. E., "Autonomous Rendezvous 
and Docking - A Commercial Approach to On-Orbit Technology Validation", Presented at the 
Autonomous Rendezvous and Capture Review, Williamsburg, VA, November 21, 1991. 
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MODELING AND EXPERIMENTAL VERIFICATION 
OF SINGLE EVENT UPSETS 


T.N. Fogarty, J.O.Attia and A. A. Kumar 
Laboratory for Radiation Studies 
Prairie View A&M University 

and 

T.S.Tang and J.S. Linder 
Texas A&l University 


Abstract 

This paper reviews the research performed and the results obtained at the 
Laboratory for Radiation Studies, Prairie View A&M University and Texas A&I 
University, on the problem of Single Events Upsets, the various schemes 
employed to limit them and the effects they have on the reliability and fault- 
tolerance at the systems level, such as robotic systems* 

Introduction 


Random access memory (RAM) based on CMOS technology has gained wide 
acceptance in space applications [ 1 ] ( 2 ] . It is known that CMOS Static RAMS 
show an upset sensitivity to single energetic heavy ions including gold, 
krypton and bromine which is called a single event upset (SEU). Immunity 
from SEU errors caused by protons or heavy ionizing particles is a requirement 
for reliable spaceborne integrated circuits. Computer simulation is an 
important mean to predict, analyze and verify the affects of SEUs on SRAMs. 
Figure la shows the circuit description of a six-transistor SRAM cell which 
has been used in SEU analyses and computer simulations. When the logic state 
is set such that node A is biased at Gnd and node B is at V*, the drain 
junctions °f Ml and M4 are sensitive regions. If an ionizing particle hits 
the junction of Ml, holes will be collected, resulting in a positive voltage 
spike at node A called the n-hit. Similarly if an ionizing particle hits the 
junction of M4, electrons will be collected, resulting in a negative voltage 
spike at node B, called the p-hit. If the voltage spikes are of sufficient 
amplitude and charge neutrality cannot be established fast enough through the 
on transistor, the flip-flop may regenerate and a bit error will occur [2]. 


Upset rate in CMOS RAMS can be reduced in two ways, either the charge 
collection capability of the memory can be degraded, or cell design can be 
altered to require greater critical charge for upset. Diminished coupling 
between the inverter pair of basic RAM cell decreases the probability of 
logic upset by slowing the feedback of gate voltage variations. If the gate 
vol age of the hit inverter remains stable during current impulse, the hit 
inverter will reestablish its prehit logic state. Maintenance of the hit 
nyerter gate voltage at or near the prehit level is accomplished by 
maintenance of the voltage at the opposite information node. Thus the upset 

k V f ty 18 dec ? eased by inc reasing the time constant of the feed back 
paths between the inverters by adding feedback resistors. 


obtained at the BNL Twin Tandem Van de Graaff Single 

lb L are , presented for 16K ' 64K and 256K Rad Hard Static 
rams ( SRAM) [3] [4]. Prior total dose radiation (1 MRAD Si of 2 MeV protons) 


° Ckn0wledge R Kohler AT&TBL, P.M.Kibule <ft V.Zajic (BNL) formerly at Hampton University and our 


♦* Partially supported by NASA-JSC-NAG-9-659, 9-331 
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produces imprinting in the Rad Hard SRAM and significantly lowers the SEU LET 
threshold. In commercial, non-rad hard devices, one would expect a reduction 
in threshold LET at much lower doses. These results are compared to SPICE 
models of SRAM with Resistive Feedback to limit Single Event Upsets (SEU). 
Passive/Active resistance networks, switched capacitor networks for limiting 
SEU are modeled and compared. The MOS transistor active resistance feedback 
method of limiting SEU has the advantage that in the absence of Cosmic Ray 
induced charge, the operation of SRAM is not degraded by high resistance of 
the active network [6-17]. 

An extension of this work considers the distribution and dependency of these 
radiation induced defects to estimate the reliability and fault- tolerance of 
more complex systems such as A/D converters and CPUs (errors in instruction 
set and in active logic) [21] [22). 

BNL Twin Tandem Van de Qraaff SEU Test 

The beam energy and LET are monitored by four silicon surface-barrier 
detectors calibrated by an Am-241 alpha source, which makes the energy 
measurement independent of the information from the accelerator control room. 
The irradiation chamber contains an adjustable iris aperture and a mobile test 
board. The effective LET can be varied by rotating the board on the vertical 
axis to change the particle angle of arrival (Figure lb). The effective LET 
is given by: 


LET 0ff = 


_ LET 


Cos0 


where theta is the angle of incidence of ions. A unique feature of the 
facility is laser optics for accurate device positioning into the beam. A 
neon laser is placed in the beam line, so that the beam spot can be simulated 
on the test board by visible light before the experiment. Evacuation and 
ventilation procedures are fully automated. Device positioning, beam 
diagnostic, data collection, and data management are computer controlled by 
user friendly menu driven software. The system specifications are summarized 
in Table 1. 


Table 1. Brookhaven SEU Test Facility 
Specifications 


Flux 

lOLlO^ions/cmVscc 

Beam Uniformity 

90-98% over 3 in. 
diameter 

Beam Aperture diameter 

0.1-1.4 in. 

Test board work area 

6 by 9 in. 

Particle angle of Arrival 

0°-73° 

Effective LET in Si 

1.4-280 MeV cm 2 /mg 

C08t 

$520/h 

$420/h for exempt users 


Radiation Hard Technology 

The devices were fabricated by AT&T-Bell Laboratories CMOS Twin-Tub IV [4] 
"1.25um H Rad-Hard technology. As the density of SRAM increases from 16K to 
64K and 256K, decrease in lateral dimensions results in lower critical charge 
due to lower gate and junction capacitances. In order to compensate for this 
increase in SEU susceptibility, the following process modifications were made; 

1) Oxide thickness was reduced and doping density was raised to increase 
specific capacitances. 

2) A twin-Tub process on thin epitaxial substrates reduced the charge 
collection by cutting off the funnel. 

3) Feedback polysilicon resistors were used to reduce the SEU rate by 
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longer decoupling times for the two inverters. 

In addition, radiation hardened gate oxide minimized threshold voltage shifts 
and transconductance degradation, radiation hardened field oxide eliminated 
parasitic leakage currents, and modified light doped drain (LDD) N-channel 
transistors minimized hot carrier effects. 

A "2um" design rule 16K, and three "1.25um" design rule (a 64K , a 256K and a 
32K x 8 ) were available for our experiments. The feedback resistor values at 
room temperature were 82, 109, 151, and 240 Kohms for the "2um" SRAM and 520 
ND 670 Kohms for "1.25um" SRAMS. Other important processing and design 
parameters are given in Table 2. 


Table 2. Processing and design parameters 



"2 um" 

"1.25 um" 

Channel width/length [urn] 

NMOS 

5.0/2.0 

1.50/1.25 

PMOS 

4.25/2.0 

1.75/1.50 

Drain area [um 2 ] 

NMOS 

50.3 

8.9;7.0 

PMOS 

25.5 

4.3;5.4 

Gate oxide thickness [A] 

215 

240 

Epitaxial thickness [um] 

1.7 

1.7 

Surface doping density [#/cm 3 ] 

N -substrate 

4x10'* 

4x10“ 

P-well 

1.5xl0 17 

1.5xl0 17 


SEU Test Experimental Set-up 

Two heavy ion beams were employed in our experiments, m Br and m Au . Most of 
our measurements were performed at elevated temperatures 80-125°C. An 
individual temperature controller described in [5] was used for each DUT. The 
chip was plugged into a home-made DIP socket with some extra pins on each 
side, and thermally coupled with a power resistor from underneath with high- 
temperature epoxy or heat-sink grease. The temperatures were measured by a 
precision ACE-48006 thermistor attached with the same compound from the top. 
The required power was less than 10W for 20 pin SRAM chips, and the nominal 
power of the heating element was even smaller than that. In the experiment, 
temperature was controlled within +1°C. 

The experiment was done using a MOSAID Memory Tester with 256K memory depth. 
Two testing modes are possible with this memory tester, static and dynamic. In 
static testing a bit map of errors is available but multiple hits are not 
recorded. In dynamic testing, the memory is checked for SEU during 
irradiation. If an error is found, the error counter is incremented by one 
and the error is corrected. Consequently, the error map is no longer 
available as it was in static testing. 16K, 64K & 256k CMOS SRAMs with 
feedback resistors performed satisfactorily in our test fixture down to 3.5 V. 
in both static and dynamic testing. However, preliminary measurements 
indicated a considerably higher SEU cross section in dynamic testing (Figure 


Because of cable capacitance (8 ft of flat cable between the memory tester, 
vacuum feedthrough, and the DUT), the constant memory reading/ loading in 
dynamic testing resulted in 1 V noise on the power supply line. The lower 
average power supply voltage decreased the critical charge and increased SEU 
sensitivity. It is also seen in Fig. 2 that dynamic testing performed at 5.0 V 
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exhibited the same critical LET as static testing at 4.5 V. By adding a fast 
tantalum capacitor 0.47 uF close to the DUT, the noise was reduced to 0.1V and 
the SEU cross section decreased. Unfortunately, we did not collect enough 
data to determine the critical LET. While dynamic testing is closer to device 
field operation, the measured SEU cross section can be substantially varied 
depending on the capacitive coupling between power supply and ground. To 
avoid this variation, all subsequent measurements were performed in the static 
mode. The power supply voltage was 5.0V and 4.5V for 16K and 64K SRAMs, 
respectively. 

Ufoftp a single memory cell is upset twice during the static testing, no error 
Im jnecorded. If the same memory cell is upset one more time, only one error 
1 M recorded and so on. Clearly, the observed number of SEU is smaller that 
tbm actual turn iber due to the possibility of multiple upsets of the same memory 
cell* D e p ending on the parity, the multiple upsets either escape observation 
entirely or are recorded as a single SEU. The following correction was 
applied to all experimental data obtained by static testing to account for 
multiple upsets: 


P*ga - & MB “ PsBO + PsEV “ • • • • 


^SBU 


where the probabilities p'^ and pja, are given by the observed and corrected 
numbers of SEU, respectively, divided by the memory size (16K or 64K) . In 
order to keep both this correction and statistical uncertainties small, it is 
good practice to make the number of SEU approximately 10% of the memory size 
by accumulating an appropriate fluence. 

SEU Test Results and Discussion 

Mo failures of any of the tested 16K SRAM (TA670) were observed at room 
temperature up to the effective LET =* 160 MeV cnr/nig with any value of the 
feedback resistors in either static or dynamic testing. The same statement 
applies for testing at temperatures up to 50°C. 

Data obtained for four 16K SRAMs with variable feedback resistors at 110°C .■> 

with Br ions using MOSAID Memory Tester are shown in Figure 3a. Data obtained 
for 64K SRAMs at the same temperature and with the same ions using the MOSAID 
Memory Tester are displayed in Figure 3b. The critical LET (defined as LET at 
which the SEU cross section drops by a factor of 500 compared to its saturated 
value, i.e. approximately 10-5 cm 2 for our 16K SRAMs) increased with the 
increasing value of the feedback resistance as expected. Since both feedback 
resistors values and the power supply voltages were different for 16K and 64K 
SRAMs, we find little sense in comparing results for the two design rules. 

All feedback resistor values are given at room temperature. They can be 
calculated at elevated temperatures using a known temperature coefficient for 
polysilicon resistors. The coefficient was measured for a resistor on the 
test chip made by the same technology (see Figure 4). The dependence was 
found to be exponential, where k is the Boltzman constant and -E = 0.0899(10) 
eV the activation energy. 

The SEU cross sections of the 16K SRAM with 82 Rohm feedback resistors were 
measured at three different temperatures, 90, 110, 125°C, and with two ion 
beams, Br and Au. Results obtained with the HP-8180/82 Data Generator/ 
Analyzer and with the MOSAID Memory Tester are shown in Figures 5a and 5b, 
respectively. A difference was observed in the critical LET for the two ions. 
The difference is, at least partially, caused by energy loss in the 
passivation layer which has an equivalent thickness of approximately 1 mg/cm 
of Silicon. The LET vs. energy curves for Br and Au ions imported from 
Ziegler's Tables are plotted in Figure 6. A small energy loss shows 260 MeV 
81 Br ions gradually approaching a maximum, thus the LET remains relatively 
constant. On the other hand, 310 HeV I97 Au ions have already passed the maximum 
and are in the region of steep decay, thus a significant LET reduction is 
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expected. For example, in our second experiment using the MOSAID Memory 
Tester, the critical LET measured at 110°C was found to be 63 and 80 MeV cm 2 /mg 
for 260 MeV 8I Br and 310 I97 Au ions respectively, thus reducing the difference 
to 10 MeV cm 2 /mg. Preliminary measurements showed that LET in Ziegler’s 
Tables is underestimated by 8% for Br and by 4% for Au in the energy regions 
of interest. Such a correction would reduce the critical LET difference to 8 
MeV cm 2 /mg. The most probable explanation of the residual difference is based 
on charge collection effects, where the ion track of Br in the device 
sensitive region is longer and more diffuse than that for Au, thus there is 
less recombination and more charge is collected. 

Using the error map, the SEU cross sections measured at high angles of 
particle arrival were corrected for the chip package shadow (up to 40% and 6% 
at 0-73° for 16K and 64K, respectively) but they still tapered off. The most 
drastic effect was observed for the 16K SRAM with 240 Kohm feedback resistors 
tested at 125°C with Au ions (see Figure 7). A sharp maximum occurred around 
the effective LET « 143 MeV cm 2 /mg and the SEU cross section dropped to zero 
above 180 MeV cm 2 /mg. The corresponding angle of incidence were 56° and 63°. 
After penetrating the passivation layer, the effective LET became 121 and 139 
MeV cirr/mg at these two angles. The calculation proves that the energy loss 
itself is insufficient for the explanation of the observed data. At high 
angles of incidence, the collected charge is probably shared by two or more 
neighboring nodes and becomes insufficient to upset either one. 

After initial SEU testing, one of the 64K SRAMs was exposed to 2 MeV protons. 

A total dose of 1.3 MRad(Si) was accumulated, 0.65 MRad(Si) without operating 
bias and 0.65 MRad(Si) in the memory state "all 0". After proton irradiation, 
the SEU test was repeated in both "all 0" and "all 1" memory states. While 
SEU cross section in the state "all 1" did not differ from its pre-radiation 
value, the SEU cross section in the state "all 0" showed a slight increase 
(see Figure 8). The SRAM was found to prefer the state in which it was 
irradiated. Ionising radiation induces bias dependent threshold voltage 
shifts and nobility degradation which cause a CMOS SRAM cell imbalance. Since 
the most sensitive strike location, for the present technology, is the OFF P- 
channel drain which is restored through an N-channel transistor, the cell 
imbalance is defined as the difference between N-channel threshold voltage 
shifts. As the 64K, 256K and 32K x 8 SRAMs utilize the same "1.25um" process 
and the same device geometry one would expect similar SEU response. The 64K 
and 32K x 8 do in fact show similar response (see Figure 9). However, a 
greater feedback delay time for the same LET threshold is apparent in the 
256K. This can only be attributed to circuit pattern effects [4a]. 


Hardening Approaches for CMOS SRAM 


Resistive Hardening 

SPICE simulation of SRAM cells with feedback resistors between the inverters 
[12] [13] [4a] shows agreement with critical LET threshold considering the 
decrease in delay as the polysilicon resistor value decreases with temperature 
increase thus decreasing LET threshold. Cell write delay times are shown in 
Figures 10 and 11 for the "1.25um" and "2um" processes respectively. The SRAM 
cells with feedback resistors are shown in Figure 12a. For the "2um" process 
SRAM disturbed by a 1 mA exponential pulse, the results at room temperature 
and 87°C are summarized in Table 3. 


Table 3. 16K SRAM, 1mA Exponential Pulse 


Temp. 

Resistance 

Max. Wid.to Q, 

27°C 

82K 

0.8nsec 

87°C 

46K 

0.3nsec 


Integrating the current pulse one finds an estimate of the disturbing charge 
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necessary to reach critical charge and cause an SEU in the SRAM cells at room 
and elevated temperatures. Clearly less charge is needed to cause an SEU at 
elevated temperatures. 

Capacitive Hardening 

In the MOSIS "2um" process the specific resistance of the polysilicon 
resistors is much lower and simulation of the distributed RC feedback shows 
that the resistor value may be decreased for equivalent LET threshold. However 
because of area considerations most commercial vendors will use higher 
specific resistance and shorter polysilicon resistor length negating the need 
for this correction [14]. 

Simulation of the insertion of a capacitance between the drain to gate nodes 
shows equivalent LET threshold and superior speed when compared with the 
feedback resistor approach. In this approach the capacitance is not in the 
write path of the cell. The write time is increased for this approach as the 
inserted capacitance increases the node capacitance. But this increased time 
is much smaller than the resistive hardening concept (see Figure 12b) [15]. 

The critical charge is 4.8 pC for a capacitance of O.lpF for a rectangular 
current pulse of 3ns width and 1.6mA amplitude. For a 100 Angstrom oxide 
thickness the area for this capacitance is nearly equal to the area of an NMOS 
transistor. Speed performance becomes a problem as larger scales of 
integration are required. We conclude that the increased area for radiation- 
hardening can be sacrificed for better speed performance. 

Active Hardening 

A CMOS Transmission Gate (TG) exhibits a nonlinear current-voltage 
characteristic when it conducts; hence called nonlinear active resistor. The 
TG resistance strongly depends on its terminal voltage: increasing rapidly as 
the terminal voltage increases. This phenomenon can be utilized to increase 
SEU immunity of a SRAM cell. A SEU hardened CMOS SRAM cell using TGs as 
feedback resistors is shown in Figure 12c. In this cell, the inverter pair is 
decoupled by two TGs whose p-channel and n-channel transistors are 
respectively gated by the ground and power source. These transistors provide 
the resistance needed for increasing critical charge of the cell and also 
introduce additional capacitance to the sensitive nodes and feedback paths of 
the cell, which can effectively increase SEU immunity of the cell [16]. 

Operation of the cell can be described briefly as following. When the cell 
operates normally, the resistance of the two TGs is very low since the 
voltages across the gate terminals are very small. The cell is essentially an 
unhardened one. When one of the sensitive nodes is hit by an ionizing 
particle, electrical charges are collected at the hit node, causing a sudden 
voltage increase or decrease at the hit node while the voltages at other 
nodes are relatively unaffected. In response to the voltage increase across 
the terminals of the TG connected to the hit node, the resistance of the TG 
becomes very high. The high feedback resistance protects the stored cell data 
from SEU. 


Effectiveness of the new SEU-hardening technique was studied numerically. The 
current induced by a particle hit was simulated by an exponential pulse. The 
rise and fall time constants of the exponential pulse were set equal to 0.01ns 
and 0.25ns respectively. 

Simulations showed that the new technique improved the SEU immunity of a SRAM 
cell effectively. With TGs whose channel length was 12um and width was 1.2um, 
a CMOS SRAM cell did not upset from a current pulse with amplitude of 10mA 
and width of 1.4ns, as seen in Figure 13a. An estimate for this TG channel 
resistance is about 100 Kohm. Prevention of upset from the same current pulse 
required a passive feedback resistor of 205 Kohm. Figure 13b shows the 
simulation result. We believed that it was the combination of distributed 
channel resistance and parasitic capacitance, mainly gate capacitance, making 
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TGs more effective than lumped resistance in increasing critical charge of the 
cell. This observation agreed with our findings in [14J where distributed RC 
feedback was shown to be better than lumped RC and pure resistive feedback. 

Switched Capacitor SRAM 

It can be shown that if clock frequency is high enough, the combination of 
switches and capacitor can replace a resistor that is dependent only the clock 
frequency and capacitor [18] [19] [20] . 

The switched capacitor network employed in this work is single-phase grounded 
switched capacitor shown in Figure 14a. it consists of two switches and a 
capacitor. For SPICE simulations, the equivalent circuit that can be used for 
the single phase switched capacitor network is shown in Figure 14b. 

The switched capacitor network was implemented using MOS technology. The 
circuit is shown in Figure 15. PSPICE simulations were performed. Samples of 
the output are shown in Figures 16,17 and 18. 

The rise time, fall time, propagation delay, and time shift (T-shift) [i.e. 
the time required for the output to occur when the input signal is applied]. 
Table 1 shows the switching times for SRAM with no feedback resistors, 
switched capacitor SRAM and SRAM with feedback resistors of 40K, 80K. and 

i c cw 9 9 


From Table 4, it can be seen that the rise time, fall time, propagation delay 
and T— shift of the switched SRAM and SRAM without feedback resistance are 
small compared to those of SRAM with feedback resistors. In addition, 
switching times of the switched capacitor SRAM are similar to those of SRAM 
without feedback resistors. 


Table 4. Switching Times of SRAMs 


reedback 

Pall Time 

Rise Time 

Prop, dly 

T-Shift 

Value 

* E-9 

* E-9 

* E-9 

* E-9 

SRAM w/out 

feedback 

0.4087 

0.380 

0.395 

0.765 

SC SRAM 

0.280 

1.029 

0.654 

0.769 

40K 

1.70 

1.581 

1.64 

2.296 

80K 

3.15 

3.121 

3.135 

4.019 

150K 

5.73 

5.503 

5.616 

6.972 


Systems Level Analysis 

This section addresses the problem of reliability and fault-tolerance from a 
systems point of view. Due to space considerations we will only provide a 
road description of our approach and model. We view the system in a 
hierarchical fashion, viz., the system is considered as a collection of 
unctional units, each functional unit made up of several functional subunits, 
and so on to the primitive elements (Figure 19). The black circles in the 
f, Sf a , 5*5*5 "if reduc ible modules" — subsystems that cannot be further 
modules " d int ° The white circles refer to "reducible 

! ny w 9 f Ven description, the system is seen to have a set of modules 

*¥*?* ” or instance, a private branch exchange system (PBX) consists 
_5 ■•quential" and "parallel" subsystems (see Figure 20). Another 

xampie is the Intel 80386 microprocessor consists of nine logical units t bus 
“ nt J r , ca ' prefetch unit, instruction decode unit, execution unit, 

unit, protection test unit, segmentation unit and paging 
unit. The execution unit is further divided into ALU and the 32-bit register 
tiles. Each of these in turn can be broken down into gates and then to 
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transistors. Two important aspects of our approach are the inclusion of: 

1. "correlational dependence' ' of failure rates; for instance, the 
occurrence of one failure may accelerate the occurrence of another; and 

2. explicit modeling of failure rates in terms of the geometrical, circuit 
and material parameters. 

Each such unit is then modeled in terms of rate equations that would allow a 
computation of the MTTF. The reliability at any given level will be 
calculated by a generalized Markov model. The knowledge of MTTF for each 
subsystem may then be used in two ways: (i) from the parametric dependence of 
the MTTF, to choose different material combinations to maximize the MTTF; and 
(ii) to replicate those subsystems with a lower MTTF, so that the entire 
system may degrade gracefully. This approach is obviously better than 
assigning an average MTTF for the whole system. 

Formulation 

The general approach in terms of a non-Markovian model is given elsewhere [21- 
23]. In the following we discuss a simplified version, the so-called memory- 
less or Markovian approximation. In this limit, which is the most common 
limit used in the literature [24-25]. The corresponding equations are: 

dP x 

— m ~ a u+i p i +fli+u p t*i +a i-u*i-i -Pv-iPi -n,/Pi. 


for i * 1, 2 , • • • ,N-1, and 


dP 0 dPu « 

at tm ® i»0 

where P f is the probability that the subsystem is in the state S t , a^ is the 
transition (failure) rate from state S { to state S jf by is the transition 
(repair) rate from state Sjto state S* and g^ f is the transition (failure) rate 
from state Si to the fatal state S f . S f denotes the "fatal failure state." 

Solution 

The above equations may be solved iteratively to yield the reliability R(t) 
and the MTTFi 


«(<)= 2^.(0 

i*0 

s(t )=c'[J!W]=c'[i , P i w] 


MTTF = S Pi(s) 


(s is the LaPlace Transform variable). 

Sample Application 

Let us apply this general solution to the case of a system with triple- 
modular- redundancy (TMR, see Fig. 21). That is, N=3. In this case, units 
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A, B and c are identical units, while unit 0 may be an identical or a 

unit. The system will, function so long one of the former three 
»w it i ff* fun f t «°” in but th * mom * nt unit 0 fails, the whole system fails. 
The failure of D is what is termed fatal or common-cause failure, ran 
electrical realization this system is a typical buffer between two large 
circuits. A, B and C could be simple inverting buffers (INRBs) while D may be 
a super-inverting buffer (INRBS).] The MTTF for this case is given by, Y 


ifrrFm A t A 2 u -t-eo.iAi-t-aa iai , 

AoAjAj -AjOojA.0 


Special Cases 



MTTF ““ 

Comment ^ 

WB3&B1 

3/a 

Three times the survival rate of a single component (or 
module), as should be expected 

a = y 

7/8a 


«0,1 * a - a l,2 “W*.t*2 > 3*K 2 a 
T»,4 »a,Vi 


k=i, mi ir=//aa, previous result 

k> 1, the MTTF is much smaller. k=2, MTTF=ll/15a 


Failure Rate Models 


Failure Phenomenon 

Failure Rate 

Electromigration 

o-AJ“«p(-E t /k B T) 

Time-Dependent Dielectric Breakdown 

| «»AQjjE|e*p(-B/E e ) 

Thermal Breakdown 

P f 

(sKpCpj^^Aj^-To) 


Similar rates may be constructed for single-event-upsets, atuck-at-faults. 
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Summary and Conclusions 


1* No SEU of AT&T RAD-HARD SRAM was observed at Room Temperature. 

2. Critical LET decreased with increasing temperature and/or decreasing value 
of feedback resistors. 

3. Critical LET was different for Bromine and Gold ions because of Zeigler 
Curve effects and the heavy ion track structure. 

4. Imprint of the memory pattern after 1.3 MRad SI TID of protons was consistent 
with the threshold voltage shift of NMOS transistors. This decrease in 
critical LET threshold is expected to be significant for non-RAD-Hard 
devices at much lower prior total dose radiation. 

5. SPICE simulation of SRAM cells with feedback resistors between the 
inverters shows agreement with critical LET threshold considering the 
decrease in delay as the polysllicon resistor value decreases with 
temperature increase thus decreasing LET threshold. 

6. In the MOSIS "2um" process the specific resistance of the polysilicon 
resistors is much lower and simulation of the distributed RC feedback shows 
that the resistor value may be decreased for equivalent LET threshold. 
However because of area considerations most commercial vendors will use 
higher specific resistance and shorter polysilicon resistor length negating 
the need for this correction. 

7. Simulation of the insertion of a capacitance between the drain to gate 
nodes shows equivalent LET threshold and superior speed when compared with 
the feedback resistor approach. However this approach is not area 
conservative. 

8. The active resistor(TG) SEU-hardening technique for CMOS SRAMs has been 
identified. This technique is effective in improving SEU immunity and needs 
no modifications of the fabrication process. The new technique shows low 
resistance except when necessary to limit cosmic ray induced charge. 

9. The switched capacitor SRAM, implemented using MOS technology, has 
characteristics similar to those of CMOS SRAM without feedback resistors. 
The switching times of the switched capacitor SRAM are comparable to those 
of SRAM without feedback resistors. In addition, it was found that the 
switching times of the switched capacitor SRAM are superior to those of 
SRAM with feedback resistors. This work shows that switched capacitor SRAM 
is a viable alternative to SRAM with feedback resistors for SEU immunity. 

10. An attempt to formulate a unified framework to compute the reliability of 
a system in the presence of fatal faults and redundant elements. In 
particular it is showed that, the inclusion of redundancy does not 
necessarily enhance the reliability of the system. Though the discussion 
has been in the framework of electronic systems, our formulation may also 
be used to describe distributed and parallel processing systems. In such 
cases, "failure" may be interpreted as "non-availability" of a processor 
(perhaps owing to its being accessed during its computation cycle) and MTTF 
might be an estimate of the length of computational time required for a 
given computation • 
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Figure 2. Effect of power supply 
voltage on SEU of 16K SRAM with 82 
Rohm feedback resistor (90*C, Br 
ions ) . Zero SEO are displayed on the 
LET axis. 
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Figure 3 a. 3X0 of 16K SRAMs with 

variable feedback resistors <110*C, 
Br ions ) • Data from MOSAID Memory 

Taster . 
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Figure 3b. SEO of 64K SRAMs with 
variable feedback resistors (110%, 
Br ions). Data from MOSAID Memory 
Tester. 
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Rohm feedback resistors at elevated 
temperatures for Br and Au ions. 
Data from HP Data/Generator /Analyzer. 
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Figure 5b. SEU of 16K SRAM with 82 
Rohm feedback resistors at elevated 
temperatures for Br and Au ions. 
Data from MOSAXD Memory Tester. 
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Figure 6. LET vs energy for "Br and 
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Figure 8. SSU of 64K SRAM with 520 
Kohm feedback resistors before and 
after proton irradiation (80°C, Br 
ions ) . 
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Figure 10. Feedback delay for 32K x 8 
("1.25um" process) cell produced by 
Cell Feedback Resistance as 
determined from write time method. 
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Figure 14a. Single phase grounded 
switched capacitor. 



Pigure 13b. Simulation result of SRAM 
with passive feedback resistors. 



Figure 14b. Continuous- time domain 
equivalent circuit with 
implementation of lossless 
transmission line. 
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Figure 16. Transient analysis of CMOS 
SRAM without feedback resistor. V(l) 
is input and V(3) is output. 


Figure 17. Transient analysis of 
switched capacitor SRAM. V(3) is 
input and V(5) is output. 



Figure 18. Transient analysis of CMOS 
SRAM with feedback resistor of 80 
Kohms. V(3) is input and V(S) output. 



Figure 19. Hierarchical 
v ^-* w of a system. 


Figure 20. A private branch 
exchange system (PBX). 


Figure 21 TMR System. 
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Summary and Conclusions 

We describe the use of simulation and contrast it to 
analytical solution techniques for evaluation of analytical 
reliability models. We also discuss the role importance 
sampling plays in simulation of models of this type. We 
next describe the simulator tool we use for our analy- 
sis. Finally, we demonstrate the use of the simulator tool 
by applying it to evaluate the reliability of a fault toler- 
ant hypercube multiprocessor intended for spacecraft de- 
signed for long duration missions. We use the reliability 
analysis to highlight the advantages and disadvantages of- 
fered by simulation over analytical solution of Markovian 
and non-Mar kovian reliability models. 

1. INTRODUCTION 

Recent work in the development of reliability analysis 
tools has produced a number of software packages that al- 
low complex system behavior to be expressed with analyt- 
ical models. The systems to which these modeling meth- 
ods are applied often are complex fault tolerant comput- 
ing systems designed for very high reliability. However, 
these systems can exhibit certain types of system behav- 
ior that require analytical models for which feasible ana- 
lytical (numerical) solution techniques are not currently 
available. In these situations the existing analytical mod- 
eling framework may be enhanced to allow simulation of 
the analytical model (i.e. a fault tree or Markov model) 
as a replacement solution method to the traditional an- 
alytical solution techniques for the model. This is the 
approach that we follow in this paper. 

The very large number of trials needed to obtain sta- 
tistically significant results historically has been a signifi- 
cant problem for the use of simulation to model complex, 
highly reliable fault tolerant systems. Recent efforts to 
overcome this problem have produced new modeling tools 
capable of obtaining acceptable results with a reasonable 
number of trials through the use of a variance reduc- 
tion technique called importance sampling. New modeling 
took which incorporate this technique have been designed 


to be compatible with the Hybrid Automated Reliabil- 
ity Predictor (HARP) modeling tool[9], which is itself a 
component of the HiRel package of reliability modeling 
tools[l]. HARP solves the same types of models as the 
simulator, but uses analytical (numerical) solution tech- 
niques instead of simulation. 

As is often the case, the development of the new mod- 
eling tool we describe here was driven by the needs of a 
specific reliability analysis project: the use of hypercube 
multiprocessors for highly reliable guidance, navigation, 
and control (G,N,& C) systems for long duration manned 
spacecraft. We are interested in exploring the use of a 
fault tolerant hypercube architecture that can use either 
hot or cold spares. It is clear from preliminary studies 
that the use of hot and cold spares with the traditional 
constant failure rate model will not meet the high relia- 
bility requirement for long duration space missions with- 
out onboard repair[ll, 12, 19]. Recently acquired empiri- 
cal data provide convincing evidence that decreasing fail- 
ure rates are common in spacecraft applications[10]. For 
these reasons, we want to be able to include decreasing 
failure rates in our reliability analysis. The inclusion of 
decreasing failure rates with cold spares requires the use 
of a non-Markovian reliability model which is substan- 
tially more difficult to solve analytically than a Marko- 
vian model that assumes constant failure rates. Given 
the current state of the art, analytical solution of such 
non-Markovian models generally is tractable only for very 
small simple models, whereas the model of the above hy- 
percube system is very large. The cumulative effect of all 
of these factors led us to the use of simulation modeling. 

In this paper we summarise the use of simulation as a 
modeling method and describe how it can be applied to 
the evaluation of analytical system models. We compare 
evaluation of'analytical models by simulation to evalu- 
ation by analytical solution techniques and describe the 
role of importance sampling in our implementation of sim- 
ulation. We next describe the simulator itself and the 
process of specifying a model for use with it. We then 
illustrate the use of the simulator by applying it to a hy- 
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percube architecture proposed for a G,N,& C system for 
long duration spacecraft. We explore the effect of assum- 
ing decreasing failure rates for active and cold processors 
within the hypercube instead of constant future rates, 
and demonstrate the advantages that simulation provides 
over analytical solution methods for such system models. 

2. SIMULATION MODELING FOR RELIABILITY 
PREDICTION 

The usual method of using simulation to evaluate re- 
liability and performance of systems involves building a 
computer model of the system, generating events of in- 
terest (i.e. component failures), and observing the re- 
sponse of the model to the generated events. The timing 
and types of events are generated using probability dis- 
tributions which are assumed to govern event occurrence. 
Values are sampled from the appropriate probability dis- 
tributions and are used to specify which type of event oc- 
curs next and when that occurrence will be. A sequence 
of events is generated in this manner until either the mis- 
sion time expires or the system fails. Such a sequence of 
events provides one instance of how the system would be 
expected to behave in the environment characterized by 
the governing probability distributions and is referred to 
as a “history” or “trial”. The model is evaluated at the 
end of a trial to determine measures of interest such as 
whether the system is still operating (reliability) or how 
much work was accomplished (performance), etc. This 
process is then repeated numerous times to obtain aver- 
age values for the measures of interest and accompanying 
sample standard deviations. From probability theory it 
follows that as the number of trials increases, the average 
value obtained in the simulation approaches more closely 
the actual value that characterizes the long run behavior 
of the system as expressed by the model. The standard 
deviation, which is a measure of the expected closeness of 
the simulation average to the actual value, is proportional 
to (where n is the number of trials)[18]. Hence ob- 
taining a highly accurate value for a measure of interest 
may require a very large number of trials. 

2.1 Analytic Solution Methods vs. Simulation 

An alternate approach to reliability evaluation involves 
building an analytic (mathematical) model to express the 
relevant behavior of the system. A number of different 
analytical model types are in widespread use. One very 
successful analytical model type is the Markov chain and 
its generalizations (non-Markovian discrete state models). 
These models express system behavior by identifying a 
number of distinct states in which the system may be. 
The system can be in only one state at a time, and from 
time to time makes a transition from one state to an- 
other. The distribution of the time the system spends 


in individual states and the characteristics of the tran- 
sition rates between states differentiate Markovian and 
non-Markovian modeis[13]. Analytical models are usually 
solved using either direct or numerical methods, so often 
they can give answers with greater accuracy than simular 
tion methods for a comparable amount of computational 
effort. However, analytical solution methods suffer from 
requiring much more memory storage for data structures 
than simulation methods. As a result, models that be- 
come too large to be accommodated by analytical solution 
methods might still be within reach of simulation tech- 
niques. In addition, increasing behavioral complexity in 
analytical models requires analytical solution techniques 
with increasing computational requirements. Hence to 
solve a model of sufficient complexity, an analytical solu- 
tion method could require more (rather than less) execu- 
tion time than a simulation method far a comparable level 
of accuracy in the output. In cases fike these where the 
limitations of analytical solution methods are exceeded,, 
simulation provides a useful alternate approach. 

The drawback to building a computer simulation model 
of a system under study is that constructing the model 
and validating it is often a complex, time consuming, 
and error-prone process. An alternative is to apply sim- 
ulation not to a model of the system itself, but to an 
analytical model of the system such as a Markovian or 
non-Markovian model. With this approach there is of 
course the problem of constructing the analytical model. 
However, this tends to be easier than constructing a 
system-level simulation model. Also, the topic of ana- 
lytical model construction has been addressed by a num- 
ber of researchers in the past several years and tools 
have been created to assist in model construction (see [3] 
and [4] for brief surveys of tools for automated Markov 
model construction). The approach we have chosen for 
the current study applies simulation to Markovian and 
non-Markovian models of the hypercube multiprocessor 
system. This allows us to capitalize on previous work 
performed by the authors on the hypercube system using 
Markovian models[5] and permits us to extend the scope 
of that work. 

2.2 Simulation for Evaluation of Markovian and non- 
Markovian Models 

Markovian and non-Markovian discrete-state models 
can be evaluated by simulation in the following way. Each 
trial represents a single traversal path among the states 
of the model. The common beginning point for all trials 
is at an initial state in which all system components are 
assumed to b'e operating correctly. Upon entry into each 
state, the process is begun for determining the time of 
transition out of the current state and which state the 
system goes to next. The time to next transition is sam- 
pled from a probability distribution that depends upon 
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Figure 1: Hypercube Multiprocessor System 


the failure rates of the components still active. If the fail- 
ure rates of all components are constant, the model is a 
Markovian model. If the component failure rates are all 
functions of mission time (i.e. non-constant), the model 
is a non-homogeneous Markov model. If the component 
failure rates are individually functions of more than one 
time variable (i.e. there is more than one “clock” in the 
system upon which component failure rates may depend), 
the model is a non- Markovian model. We use all three 
types of models in the present study. Once the time to 
next transition has been determined, a sampling from a 
second distribution is done to determine which of the re- 
maining operating components will experience the fail- 
ure that is the cause of the transition out of the state. 
The determination of the sampling distributions is de- 
scribed in [14] and [16]. We note that this formulation 
of the simulation process can accommodate the use of 
Fault/Error Handling Models (FEHMs) to implement be- 
havioral decomposition for incorporating imperfect fault 
coverage as is done in HARP[15]. Although that capa- 
bility was available, we did not consider imperfect fault 
coverage in the present study. During each trial succes- 
sive inter-state transitions are generated until either the 
mission time is exceeded or the system fails, causing the 
trial to end. The system unreliability is then estimated 
from the proportion of trials during which the system 
failed before the mission time was reached. 

2.8 Importance Sampling 

A major characteristic of highly reliable systems is that 


system failure events are extremely rare. This means that 
a large majority of the trials are likely to end by the mis- 
sion time expiring rather than through a system failure. 
Since system failures are the events of interest, a very 
large total number of trials must be run before a sufficient 
number of system failures occur to provide a meaningful 
estimate from the proportion of failure trials to total trials 
(i.e. an estimate of the system unreliability). A variance 
reduction technique called importance sampling may be 
employed to reduce the total number of trials required. 
An excellent introduction to importance sampling may 
be found in [6]. The basic idea behind importance sam- 
pling is to select an alternate distribution from which to 
sample which has much higher probability density than 
the original distribution in the regions of interest where 
the original distribution’s density is very small. Parity to 
sampling from the original distribution is maintained by 
weighting the observations sampled from the new distri- 
bution to reflect the relative difference in density magni- 
tude between the two distributions. For example, if the 
density of the new distribution is four times greater than 
the density of the original distribution in a certain region, 
then a failure event observed in that region by sampling 
from the new distribution is counted as only \ of a failure. 
The importance sampling techniques implemented in the 
simulator we used for this study, called forced transitions 
and failure biasing, are described in [14]. Both have the 
effect of emphasizing component failure events in order to 
increase the number of trial terminations due to system 
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Figure 2: Modeling Functional Dependencies due to Pro- 
cessing Node Interconnections 



Figure 3: Fault tree model of Architecture 1 Processing 
Node with Cold Spares 

failure, hence reducing the total number of trials needed 
in order to accumulate a sufficient number of system fail- 
ure terminations to provide an acceptable estimate of the 
system unreliability. 

2.4 Simulator Description 

The original version of the simulator w£ used for our 
analysis was designed by Lewis[17] and implemented at 
Northwestern University. It required a system model to 
be described as a set of components arranged in groups. 
Each group could optionally have cold spares, and could 
have either a constant or a Weibull increasing failure rate. 
Each group could also have a Fault/Error Handling Model 
(FEHM) associated with it to allow the use of behavioral 
decomposition as is done in HARP. System failure cri- 
teria were specified in the form of a set of component 
cut sets which the analyst had to derive from a combi- 
natorial model of the system (for example, a fault tree). 
For our study, we modified this simulator to enable it to 


use decreasing as well as increasing Weibull failure rates, 
and to allow it to accept the input model in the form 
of a dynamic fault tree (see below) rather than as a set 
of component cut sets. The resulting simulator program 
accepts its input model in the same form as the HARP 
program, and accepts input files with the same format as 
HARP. In addition, it is capable of evaluating all models 
that HARP is capable of evaluating, making it completely 
compatible with HARP. This is an important advantage 
because it allows the reliability analyst to develop his/her 
system model once and then input it to whatever evalua- 
tion program is most appropriate depending on the char- 
acteristics of the model and the programs. It also allows 
a comparative evaluation of the performance of the two 
programs by applying them both to the same model(s). 

8. SYSTEM MODEL 

The hypercube multiprocessor system and the mode 
of it that we use in this study are described in [3}, and 
[5] under the name of Architecture 1. We give a brief 
description of it here. The architecture is shown in figure 
1. It consists of a 3-dimensional hypercube configured ' 
as two fault- tolerant 2-dimensional modules, each witlr a 
spare processing node. The processing nodes themselves 
are multiprocessors containing four . active processors and ; 
a spare processor. The spare processor can, be either a 
hot or cold spare. The structure of the proce3»fe^ nodes v , 
is also shown in figure 1. Each processing ^ode^onunpm^ - c 
cates with other processing nodes in the system thljcTUgh 1 
four ports. For the system to be operational all eight pro^ r;7 
cessing nodes must be operational and must all be able * s 
to communicate with each other. Therefore, the system v . 
will be considered failed if any processing node fails and 
a spare processing node is unable to take over or if any 
two nodes in the hypercube are unable to communicate 
with each other. 

Although the form of the analytical model that is ac- 
tually evaluated is a Markovian/non-Markovian discrete- 
state model, it is specified by the reliability analyst in the 
form of a dynamic fault tree[ 3, 8]. When simulation is not 
used for model evaluation, the dynamic fault tree can be 
converted into a Markov chain which can then be solved 
numerically for state probabilities. When simulation is 
used for model evaluation, the discrete-state structure of 
the underlying Markovian model is inherent in the sim- 
ulation process and the dynamic fault tree is used only 
to determine whether a state which has been entered is a 
failure state/ 

A dynamic fault tree is a generalized fault tree model in 
which the traditional set of combinatorial fault tree gates 
is extended to include several non-standard gates that 
are designed to express sequence dependent behavior. Se- 
quence dependent behavior is behavior that depends in 
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Node fails 



Figure 4: Fault tree model of Architecture 1 Processing 
Node with Hot Spares 



Figure 5: Effect of Weibull DFRs on System Unreliability 
(Hot Spares) 


Component 

Initial constant failure rate 

Shared Memory 
Intrarnode bus 
Processor 

3.477 x 10“ 7 
1.147 x 10- 7 
1.990 x 10“* 


Table 1: Initial Constant Hasard Rates (failures/hour) 
for Components in Processing Nodes 

some way on the order in which events occur. The hyper- 
cube system under study exhibits two instances of this 
type of behavior: functional dependencies (the failure of 
one component causes one or more other components to 
either fail or become unavailable) and cold spares (a cold 
spare cannot fail while it is “cold"; it can fail only af- 
ter it has been activated to substitute for a failed active 
component). The functional dependencies appear in the 
interconnections between the processing nodes; specifi- 
cally, if either the internode link or one of the two ports 
on either side of the internode link fails, the remaining 
two components (link and/or port(s)) become useless to 
the remaining operation of the system and hence may 
be considered to be effectively failed themselves. These 
functional dependencies are modeled with functional de- 
pendency gates , as shown in figure 2. Cold spares are used 
within the processing nodes and are modeled using a cold 
spare gate % an example of which appears in figure 3. 

Figure 4 models a processing node when the spare pro- 
cessor is hot (i.e. active and running from mission start 
just like the four initially active processors). The 2-out- 
of-4 gate for which the four ports are inputs reflects the 
effect of the message routing protocol[5]. Figure 3 models 
the processing node when the spare processor is cold. Di- 
agrams of fault trees modeling the full architecture were 
omitted from this paper due to lack of space. The inter- 
ested reader may find them in [3], 

l ANALYSIS RESULTS 

We evaluated the system model for the cases where 
all components had constant failure rates with hot or 


cold spare processors (time homogeneous Markov mod- 
els), various components had Weibull DFRs with hot 
spare processors (non-homogeneous Markov model), and 
various components had Weibull DFRa with cold spare 
processors (non-Markovian model). For this paper our 
primary purpose is to illustrate the use of simulation to 
evaluate the models and contrast it with analytical solu- 
tion techniques. Therefore we will use here only selected 
results from our analysis to compare the advantages and 
disadvantages of simulation vs. analytical solution meth- 
ods. A more complete reliability analysis of the hyper- 
cube system is found in [2]. Our primary analysis goal 
was to determine whether assuming Weibull decreasing 
failure rates (DFRs) for components instead of constant 
failure rates would result in a sufficient improvement in 
predicted system reliability to conclude that the archi- 
tecture was adequate to successfully complete a 10 year 
mission. Results using constant failure rates[2] indicated 
that the proposed architecture would be inadequate, with 
the probability of system failure exceeding 60% after 10 
years. Initial attempts to evaluate the model with HARP 
(which uses analytical solution techniques) were not suc- 
cessful due to the large size of the model. The dynamic 
fault tree model of the system contains 70 basic events 
(110 components total), and 175 fault tree nodes (basic 
events + gates). It produces a Markov model with many 
thousands of states. Furthermore, when Weibull DFRs 
are assumed together with cold spares, the size and com- 
plexity of the resulting non-Markovian model is well be- 
yond the capability of any analytical solver tool that ex- 
ists today, both in terms of memory and execution time 
required for its solution. In contrast, our simulator was 
able to evaluate the model with none of the problems ex- 
perienced by HARP. Components with decreasing failure 
rates were assumed to have an initial failure rate A« X p 
given in table 1 which declines monotonically over the 
mission time according to the Weibull failure rate expres- 
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Mission 

Time 

(Years) 

All Components 
Constant FRs 

Processors 
Weibull DFRs 

Processors 
and Ports 
Weibull DFRs 

All Components 
Weibull DFRs 

1 

.249 ±.016 

.0250 ± .0031 

.000519 ± .00022 

.000255 ±.00013 

2 

.271 ± .016 

.0489 ±.0048 

.00147 ±.00031 

.000361 ± .00015 

3 

.312 ±.017 

.0738 ± .0065 

.00286 ± .00044 

.000439 ±.00017 

4 

.361 ±.018 

.0988 ± .0091 

.00481 ±.00078 

.000504 ±.00019 

5 

.419 ±.018 

.126 ±.014 

.00729 ± .0013 

.000550 ±.00020 

6 

.475 ±.018 

.152 ±.017 

.0102 ± .0018 

.000638 ± .00031 

7 

.530 ±.018 

.176 ±.019 

.0135 ± .0024 

.000673 ±.00033 

8 

.576 ± .017 

.202 ±.023 

.0173 ± .0037 

.000718 ±.00036 

9 

.609 ± .016 

.231 ±.031 

.0208 ± .0045 

.000766 ±.00041 

10 

.631 ± .013 

.257 ±.036 

.0257 ±.0091 

.000777 ±.00041 


Table 2: Effect of Weibull DFRs on System Unreliability (Hot Spares) 


sion: 

* ( 1 ) 

where a is the Weibull shape parameter [20] which is as- 
sumed to have the value <x = 0.5. All components not 
having DFRs were assumed to have constant failure rates 
given in table 1. Table 2 and figure 5 show the effect of as- 
suming Weibull DFRs for various subsets of components. 
The results reported in table 2 are averaged over 10 runs 
of 10000 trials per run. The effect of assuming Weibull 
DFRs for increasing numbers of the components clearly 
results in decreasing system unreliability. The result of 
assuming Weibull DFRs for all components is a difference 
of about three orders of magnitude in the system unre- 
liability (from 0.631 ± 0.013 when all components have 
constant FRs down to about 0.777 x 10~ 3 ± 0.41 x 10 -3 
when all components have Weibull DFRs). 

The above discussion illustrates the advantage that 
simulation can have over analytical techniques: simular 
tion may be able to evaluate models that are beyond the 
reach of analytical techniques both in terms of memory 
and execution time. Furthermore, if only ballpark evalu- 
ations are desired, simulation may be able to produce the 
required results relatively quickly. Figure 6 contrasts the 
reliability predictions for the hypercube with hot spares 
assuming constant failure rates and Weibull DFRs for all 
components. The results are averaged over 10 runs, with 
each run consisting of only 1000 trials requiring approx- 
imately 4 minutes or less of clock time. With only 1000 
trials per run, the standard deviations are relatively large. 
Nevertheless, the outcome of the comparison is clearly ap- 
parent. 

However, simulation does have an important disadvan- 
tage compared to analytical solution techniques. If the ac- 
curacy of the evaluation is important, then the execution 
time required by simulation to achieve the required accu- 
racy increases rapidly and can quickly become uncompet- 
itive with that required by analytical solution techniques 


(provided the model is small enough for analytical solu- 
tion techniques to be used). Table 3, which shows the 
reliability of a single processing node in the hypercube 
and the execution time required to obtain it, contrasts 
the values obtained using HARP to values obtained us- 
ing the simulator with varying numbers of trials per run. 
Increases in the accuracy of the reliability estimate, as 
measured by the decreasing size of the standard devia- 
tion, require very significant increases in the execution 
time. The table clearly shows that it is better to use 
the analytical solver than the simulator, both in terms 
of execution time and accuracy of the reliability predic- 
tion. This result holds in general, and experience has 
shown that it is usually preferable to use an analytical 
solver whenever feasible rather than a simulator to evalu- 
ate a reliability model. In particular, whenever accuracy 
in results is important we feel that the use of a simula- 
tor generally should be a last resort to be pursued after 
analytical modeling techniques have been found to be in- 
feasible. 

5 . SUMMARY 

We have described a reliability analysis study which 
was performed to determine whether assuming of Weibull 
decreasing failure rates (DFRs) for components of a fault 
tolerant hypercube would significantly improve the 10 
year system reliability estimate over that obtained as- 
suming constant failure rates. Our results show that a 
substantial improvement in system reliability does result 
from assuming Weibull DFRs, indicating that a candi- 
date architecture that would otherwise be considered in- 
adequate instead could provide acceptable reliability after 
all. We also contrasted the use of simulation and analyt- 
ical solution techniques to evaluate Markovian and non- 
Markovian reliability models. Observations made from 
our analysis indicate that analytical solution techniques 
are preferable whenever the model is small enough and 
when accuracy of the answer is important. Conversely, 
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Figure 6: Ballpark Evaluation of the Effect of Weibull DFRs on System Unreliability (Hot Spares) 


Solver 

Reliability 

Estimate 

CPU time 
required 

HARP 

Simulator, 10 s trials/run 
Simulator, 10 4 triab/run 
Simulator, 10 s triab/run 
Simulator, 10* triab/run 

0.04468 
.04374 ± .0023 
.04455 ±.00073 
.04462 ±.00023 
.04463 ±.000073 

6.4 sec 
29.2 sec 
4 min 53.8 sec 
48 min 26.9 sec 
8 hrs 0 min 2.5 sec 


Table 3: Processing Node Model Evaluation Accuracy vs. Execution Time 


simulation is preferred whenever approximate ballpark 
answers for a large model are sufficient, or when the model 
is too large or exhibits system behavior too complex to 
be accommodated by analytical solution techniques. Fi- 
nally, we have described a simulator tool for evaluating 
Markov and non-Markovian reliability models which is 
compatible with the HARP (analytical) reliability evalu- 
ation program and is part of the HiRel package of reliabil- 
ity evaluation tools. There is a great advantage to having 
analytical and simulation tools be compatible with each 
other in this way (i.e., both using the same input models 
and files, and both providing the same analysis capabil- 
ity) because it allows the reliability analyst a great deal 
of flexibility in conducting the analysis. Solution meth- 
ods may be mixed and matched and applied in the most 
appropriate way to a single system model depending on 
the type and scope of the desired results. 
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Early Lunar Rover Mission Studies 




V. P. Gillespie 

NASA Langley Research Center 
Hampton, Virginia 
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October 1992 - Dynamic Model Tested 
November 1992 - JSC put Artemis on hold 
December 1 992 - Study Completed 
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The Robotic All-Terrain Lunar Exploration 
Rover (RATLER) — Increased 
Mobility Through Simplicity 


J. Bryan Pletta 

Sandia National Laboratories 

Albuquerque, New Mexico 



Abstract 


The Robotic All-Terrain Lunar Eiploration Rover (RATLER) 
- Increased Mobility through Simplicity 

J. Bryan Pletta 
Dept 9616 

Sandia National Laboratories 
Albuquerque, NM 87185 


A new concept mobility chassis for a robotic rover is described which is inherently simple 
with few moving parts or complex linkages. The RA TLER design utilizes a four-wheel 
drive, skid steered propulsion system in conjunction with passive articulation of the dual 
body vehicle. This uniquely simple method of chassis articulation allows all four wheels 
to remain in contact with the ground even while climbing obstacles as large as 1.3 wheel 
diameters. A prototype mobility platform has been built which is approximately 1 m^ 
with 0.5 m diameter wheels and all-wheel electric drive. The theoretical mobility 
limitations are discussed and compared with the results of field trials of the prototype 
platform. The theoretical model contrasted with measured performance is then used to 
predict the expected mobility of the RATLER design on the Lunar surface. 
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The Robotic All-Terrain Lunar Exploration Rover 
Increased Mobility Through Simplicity 


V. 
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RATLER Kev Characteristics 
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RATLER Photo 
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r ATT ,F,R Optimum Wheel Base 





RATLER Maximum Steo Height 
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RATLER Pathfinder Configuration 
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Lunar Surface Characteristics 
Slopes 
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Lunar Surface Characteristics 
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Degraded Crater Model 
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Lunar Surface Characteristics 
Minimum and Maximum Significant Craters 








RATLER Functional Requirements 
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RATLER Pathfinder Field Trial Results 
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RATLER Mechanical Design 
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A Multitasking Behavioral Control System 
for the Robotic All-Terrain Lunar 
Exploration Rover (RATLER) 


Paul Klarer 

Sandia National Laboratories 
Albuquerque, New Mexico 
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A Multitasking Behavioral Control System 
for the 

Robotic All Terrain Lunar Exploration Rover (RATLER) 


P. Klarer 


Sandia National Laboratories 
Advanced Vehicle Development Department 
Robotic Vehicle Range 
Albuquerque, New Mexico 


Abstract 

An approach for a robotic control system which implements so called 'behavioral' control within a 
realtime multitasking architecture is proposed. The proposed system would attempt to 
ameliorate some of the problems noted by some researchers when implementing subsumptive or 
behavioral control systems, particularly with regard to multiple processor systems and realtime 
operations. The architecture is designed to allow synchronous operations between various 
behavior modules by taking advantage of a realtime multitasking system's intertask 
communications channels, and by implementing each behavior module and each interconnection 
node as a stand-alone task. The potential advantages of this approach over those previously 
described in the field are discussed. An implementation of the architecture is planned for a 
prototype Robotic All Terrain Lunar Exploration Rover (RATLER) currently under development, 
and is briefly described. 
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Tour Presentations 
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Automated Structural Assembly 
Laboratory (ASAL) 


Ralph Will 

NASA Langley Research Center 
Hampton, Virginia 


A single robot arm constructs a 102 member planar truss with hexagonal 
reflector panels autonomous .y. The arm travels on as X-Y carnage and the truss 
is assembled on a rotating turntable. The truss hardware and end-effectors were 
developed in-house. Technology efforts have included automated error recovery, 
machine vision guidance, end-effector microprocessor control, path planning, 
sequence planning, artificial intelligence, automated end-effector change-out, and 
software design. 
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Intravehicular Automation and Robotics 


Kelly Willshire 

NASA Langley Research Center 
Hampton, Virginia 


A full-scale mockup of Space Station Freedom’s Laboratory Module has been 
built. A robotically controlled 7-DOF arm rides on a track for mobility in the 
mockup. The purpose is to investigate how automation and robotics technology 
can improve the productivity of SSF experiments, especially when astronauts are 
absent. The facility will address protein crystal growth and furnace experiments 
in 1993 and 1994. 
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Flight Telerobotic Servicer Hydraulic 
Manipulator Testbed (FTS HMTB) and 
Vehicle Emulator System (VES) 


Wallace Harrison and Robert L. Williams 
NASA Langley Research Center 
Hampton, Virginia 


The FTS HMTB is the ground-based testbed for the FTS arm which will 
be delivered to NASA JSC in June 1993. The HMTB and its control system 
is identical to the flight arm, except that hydraulic actuation is required to lift 
representative space payloads in 1-G. 

The VES is a six-legged hydraulic Stewart platform mechanism which is used 
to study disturbance compensation for external operations of space telerobotic 
systems. The class of devices modeled is compound manipulators, such as 
SPDM attached to the end of SSRM3, or a manipulator on a free-flying satellite. 
Manipulator arms will be mounted to a load cell on top the platform. Based on 
the inertial forces, an admittance model will drive the platform to emulate the 
motion of compliant manipulator vehicles in space. 
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Intelligent Systems Research 
Laboratory (ISRL) 


Robert L. Williams and Ed Hogge 
NASA Langley Research Center 
Hampton, Virginia 


The ISRL has developed dual arm, shared control of sensor-rich telerobotic 
systems. Several representative space tasks have been demonstrated over the 
years. Simultaneous control is possible combining hand controller inputs, au- 
tomatic position commands, machine vision guidance, force control, and l«!wr 
proximity control. Two standard PUMA arms have been used for some time, and 
two red un dant 8-axis arms are recent additions. 
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